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ABSTRACT 



An orbiting detector of infrared (IR) energy may be used to detect the rocket 
plumes generated by ballistic missiles during the powered segment of their trajectory. By 
measuring angular directions of the detections over several observations, the trajectory 
properties, launch location, and impact area may be estimated using a nonlinear least- 
squares iteration procedure. Observations from two or more sensors may be combined to 
form stereoscopic lines-of-sight (LOS), increasing the accuracy of the estimation 
algorithm. The focus of this research has been to develop a computer-model of an 
estimation algorithm, and determine what parameter, or combination of parameters will 
significantly affect on the error of the tactical parameter estimation. This model is coded 
in MATLAB, and generates observation data, and produces an estimate for time, position, 
and heading at launch, at burnout, and calculates an impact time and position. The effects 
of time errors, LOS measurement errors, and satellite position errors upon the estimation 
accuracy were then determined using analytical and Monte Carlo simulation techniques. 
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1. INTRODUCTION 



The TRW-built Defense Support Program (DSP) satellites have been the 
spacebome segment of NORAD’s Tactical Warning and Attack Assessment system since 
the early 1970’s. Using infrared detectors that sense the heat from missile plumes against 
the Earth background, these orbiting sentries detect ballistic missile launches DSP also 
detects nuclear detonations in support of Nuclear Test Ban monitoring The DSP system 
provides near real-time detection information in support of DOD’s tactical warning and 
attack assessment mission and is supported by a network of fixed and mobile ground 
stations that process and disseminate information to military commanders worldwide. The 
Cold War mission of the DSP system was to detect massive intercontinental ballistic 
missile (ICBM) attacks. United States’ response to such an attack only required timely 
and unambiguous warning. Missile flight times were much longer than the time required 
to launch a retaliatory attack Precise radar tracks could be established with enough time 
to mitigate effects as much as technology allowed. 

The current combat environment demands much more from launch detection 
satellites. The present threat is from tactical ballistic missiles (TBM’s), which exhibit 
much cooler bum and shorter thmst times, and possibly more depressed trajectories than 
those exhibited by ICBM’s. TBM’s can be launched from almost anywhere within a large 
geographical area of interest, with lofted or depressed trajectories. There may be no other 
sensors to quantify impact parameters in time to employ anti-ballistic missile (ABM) 
weapons or alert potential victims within the impact zone. 

For budgetary reasons, the United States is forced to use existing launch detection 
satellites to counter the TBM threat into the beginning of the next century. [Ref 1 ] More 
rapid extraction of more precise information from these existing, technology-limited 
satellites must be accomplished in the interim of acquiring the next-generation TBM 
detection satellite system. Anti-ballistic missile (ABM) systems, such as Patriot or Aegis, 
may be employed as ABM umbrellas in a tactical area of interest, provided they receive 
precise and timely cueing from these detection systems 



This raises the question: “How accurate is the information provided by DSP?” To 
begin to answer this question, an engineering error analysis must be performed to 
determine; 

• What are the errors present in the detection system? 

• Which errors have the greatest effect upon the accuracy of DSP output 
information? 

• What are the effects of the errors on the results (individually and collectively)? 

By modeling the algorithm used by the tactical warning system to determine trajectory 
elements and predict impact zone location of a detected TBM, it is then possible to 
introduce the inherent errors separately into the model, and then study their effects upon 
the results. The relative magnitudes of their effects upon the output and their effects upon 
the results are then determined. This is done both qualitatively and numerically 
(statistically), and the results are compared. 
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II. BACKGROUND 



The DSP system consists of one or more satellites in geosynchronous orbit and 
one or more ground receiving stations. A DSP satellite consists of a bus (spacecraft) 
segment and a payload (sensor) segment. The satellite is 10 meters long, 7 meters in 
diameter, and weighs over 2300 kilograms. [Ref 2] 




Figure 1 . Diagram of DSP Satellite [From Ref 2] 

The spacecraft segment receives and transmits all commands and data, provides 
attitude, spin-rate, and station-keeping functions, and provides and distributes satellite 
power. Solar cells mounted on the spacecraft’s cylindrical body and on four solar panels 
mounted opposite the sensor generate the electrical power. The sensor collects IR energy, 
provides onboard (initial) data processing, and supplies precise orientation data. 
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The satellite is placed in a circular, equatorial, geosynchronous orbit (Titan/IUS), 
oriented so that the telescope is pointed toward the earth, and spun along its longitudinal 
axis at six rpm. This configuration is called a “yaw spinner”. The axis of rotation of the 
satellite is normal to the earth’s surface (nadir-pointing). The satellite’s spin allows the 
sensor to regularly scan the earth and distributes the thermal load uniformly. The spin also 
allows one set of attitude control components to effect two-axis earth pointing control 
(roll and pitch) by time-sharing A counter-rotating momentum wheel keeps the net spin 
momentum near zero, thereby minimizing the gyroscopic effect due to coupling between 
the spin motion and the orbit rate with minimum fuel expenditure. 



The IR telescope is tilted from the spin axis, so that the photo-electric cell (PEC) 
array covers the radius of the earth. As the satellite rotates, the entire surface of the earth 
within the FOV is scanned by the IR detector, shown in Figures 2 and 3. 




Figure 2 Photoelectric Cell Array [From Ref 2] 



4 





Figure 3 Focal Plane Array Scanning FOV [From Ref 3] 
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Detection of ER sources is accomplished with the telescope and PEC-array 
portions of the sensor. The signal electronics onboard the satellite partially discriminate 
the targets from the background and the ground station completes the task. Location of 
the IR sources is derived from orientation of the IR telescope line-of-sight relative to the 
earth’s surface. Star sensors augment the sensor data for precise determination of sensor 
pointing direction. 

The PEC-array of the ER detector is mounted with the nadir end at the center of 
rotation of the telescope (see Figure 3). This array contains over 6000 detector cells that 
are sensitive to energy in the infrared wavelengths. As the PEC array scans the FOV, a 
cell passing across an ER source will generate a voltage with an amplitude proportional to 
the signal intensity. This voltage signal is termed an ER return, and is transmitted to 
ground processing stations after amplification and background filtering. A simplified 
diagram of the data distribution network is shown in Figure 4. 




Figure 4 Data Distribution Diagram [Ref 4] 

DSP has ground stations globally positioned to receive, process, and report 
mission data to users The Satellite Control Facility tracks the satellite, transmits 
commands, receives mission and satellite status data, and forwards the data to the Data 
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Reduction Center (DRC) for further processing. The DRC processes the data, extracts 
significant returns associated with missile launches and nuclear detonations, and generates 
event messages for transmission to the users over the Ground Communications Network. 
[Ref 2] 
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HI. THE ALGORITHM 



{ Thomas Jerardi has graciously given his permission to adapt and modify his 
unpublished paper, “TALON SHIELD/ALERT State Vector Estimation", which is the 
basis for the majority of this chapter. } [Ref. 5] 

The detection of a TBM launch is the starting point for any ballistic missile de- 
fense. DSP does this by detecting the IR radiation emitted by the exhaust plume of a 
launching missile. With detections by two or more spacecraft (stereo observations), trian- 
gulation of lines of sight can be used to more accurately estimate the boost-phase trajec- 
tory, which is then used to calculate launch position, state vector (position and velocity) at 
missile engine burnout, and impact position. Real-time knowledge of the launch position 
allows targeting of the launcher. Cueing for ABM weapons systems, such as Patriot or 
Aegis, requires timely and accurate trajectory information, which can be propagated from 
knowledge of the state vector at burnout. Impact time and position data is extrapolated 
from the state vector at burnout, and may be used for warning personnel within the target 
area. Therefore, the quality (accuracy) of the trajectory estimation process is of para- 
mount importance to Tactical Ballistic Missile Defense (TBMD). Understanding the 
algorithms and equations employed in the estimation process is necessary to assess the 
quality of the estimated parameters. 

The tactical parameter estimation process is composed of several tasks: initial 
estimate of the tactical parameters, nonlinear least-squares estimation (refinement) of 
tactical parameters, burn-out time estimation, state vector generation, and impact point 
calculation. The tactical parameters are: 

• To = Time of launch 

• L =Loft 

• <|)o = Launch point geodetic latitude 

• = Launch point geodetic longitude 

• ho = Launch point height above WGS-84 ellipsoid 

• Oo = Flight trajectory azimuth (true heading) 

These quantities will be explained in more detail in the following sections. 
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Observational data are used to calculate initial estimates of tactical parameters and 
to choose the appropriate TBM profile from a database. The profile trajectory is then 
used to calculate expected observations, which are then compared to the actual observa- 
tions. The differences are used to calculate changes to the initial estimates, and a least- 
squares iteration is run until the differences between expected and actual observations are 
sufficiently small. The result is an accurate determination of the TBM’s launch and trajec- 
tory parameters. An estimate of burnout time is then made, and the calculated trajectory 
is extrapolated to generate the state vector at burnout. This allows the computation of 
impact time and position using simple ballistic trajectory equations. Each of these tasks is 
described in detail in the following sections. 

A. OBSERVATIONAL DATA 

The starting point is the observational data, which are measurements of IR radiant 
intensities (IR returns) as a function of time from each of the approximately 6000 focal 
plane elements in the DSP satellite’s sensor. Attitude information (star sensor measure- 
ments, jet-firing data, momentum wheel tachometer data) is also included in the telemetry 
stream and used to determine spacecraft attitude history. The satellite’s position is de- 
termined by the Air Force Satellite Control Network (AFSCN). A global perspective of 
the observation geometry is shown in Figure 5. 




Figure 5. Observation Geometry [From Ref 5] 
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The discrimination of TBM IR returns from earth background is a complex proc- 
ess It will be assumed to have been effectively done, and that a table of observations 
corresponding to a single, boosting TBM is available. An example is shown in Table 1 . 



Index 


Time 

(sec) 


S/C 

ID 


Intens. 


Azimuth 

(rad) 


Elevation 

(rad) 


G.H.A. 

(rad) 


Dec. 

(rad) 


Radius 

(km) 


k 


T 


S/C 


I 


P 


ft 


gha 


5 


R 


1 


129.36 


1 


14.0 


4.061854 


0.115847 


0.174532 


0 


42164.17 


2 


130.30 


2 


23.0 


2.427020 


0.094741 


1.221730 


0 


42164.17 


3 


135.44 


3 


32.0 


2.062181 


0.142310 


1.832595 


0 


42164.17 


4 


139.36 


1 


29.0 


4.062497 


0.115971 


0.174532 


0 


42164.17 


5 


140.30 


2 


29.0 


2.427264 


0.094753 


1.221730 


0 


42164.17 


6 


145.44 


3 


40.0 


2.061969 


0.142405 


1.832595 


0 


42164.17 


7 


149.36 


1 


49.0 


4.063552 


0.116147 


0.174532 


0 


42164.17 


8 


150.30 


2 


60.0 


2.427660 


0.094746 


1.221730 


0 


42164.17 


9 


155.44 


3 


65.0 


2.061752 


0.142493 


1.832595 


0 


42164.17 



Table 1. Example Observational Data 



The index, k, runs from 1 to n, the total number of observations (typically less than 
20), and is used as a subscript for the remaining symbols to denote a particular observa- 
tion. Time, Tk, is the time of observation measured from midnight of the day of the obser- 
vation, and runs between 0 and 86400 seconds. Spacecraft Identification, S/C ID, is 
simply a notation used to identify which satellite is making which observation. The in- 
tensity, Ik, is the radiant intensity of the IR return, and is used to select the type of TBM 
being detected. Azimuth angle, 3k, is the azimuth of the line of sight of the IR return 
measured clockwise from true south, and has a range of 0 to 2n radians (0° to 360°). 
Elevation angle, t|k, is the elevation of the return measured from nadir, and has values 
between 0 and 0.175 radians (0° to 10°). Satellite position is given in spherical coordi- 
nates. Greenwich hour angle, ghak, is the angle between the satellite’s nadir point and the 
prime meridian, measured positive east. Declination angle, 5k, is the angle above or below 
the equator, measured positive north. (Greenwich hour angle and declination have been 
adopted instead of sub-satellite point longitude and latitude to avoid confusion with TBM 
latitude and longitude.) Radius, Rk, is the distance from the satellite to the earth’s center, 
measured in kilometers. Geosynchronous satellites typically stay within a few degrees of 
the equator, and have a radius of approximately 42,164. 17 km. 
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B. TBM PROFILE 



A TBM profile is a description of the nominal powered flight trajectory of the 



given TBM. An example TBM profile is given in Table 2. 



Time 

(sec) 

0 


Intensity 


Altitude 

(km) 


Range 

(km) 


Time 

(sec) 


Intensity 


Altitude 

(km) 


Range 

(km) 


36.0 


0.000 


0.000 


32 


60.6 


7.023 


3.195 


1 


36.3 


0.006 


0.000 


33 


62.4 


7.469 


3.491 


2 


36.6 


0.026 


0.000 


34 


64.2 


7.928 


3.803 


3 


36.9 


0.058 


0.000 


35 


66.0 


8.402 


4.132 


4 


37.2 


0.103 


0.000 


36 


68.4 


8.890 


4.479 


5 


37.5 


0.163 


0.001 


37 


70.8 


9.393 


4.844 


6 


37.8 


0.235 


0.004 


38 


73.2 


9.911 


5.229 


7 


38.1 


0.322 


0.010 


39 


75.6 


10.444 


5.633 


8 


38.4 


0.423 


0.020 


40 


78.0 


10.992 


6.057 


9 


38.7 


0.537 


0.036 


41 


81.2 


11.556 


6.502 


10 


39.0 


0.666 


0.058 


42 


84.4 


12.136 


6.969 


11 


39.5 


0.809 


0.087 


43 


87.6 


12.732 


7.459 


12 


40.0 


0.965 


0.124 


44 


90.8 


13.345 


7.973 


13 


40.5 


1.136 


0.171 


45 


94.0 


13.975 


8.511 


14 


41.0 


1.321 


0.226 


46 


96.0 


14.622 


9.075 


15 


41.5 


1.520 


0.292 


47 


98.0 


15.288 


9.665 


16 


42.0 


1.733 


0.367 


48 


100.0 


15.972 


10.282 


17 


42.5 


1.962 


0.453 


49 


102.0 


16.675 


10.928 


18 


43.0 


2.204 


0.550 


50 


104.0-^- 


17.397 


11.604 


19 


43.5 


2.460 


0.658 


51 


104.6 


18.140 


12.309 


20 


44.0 


2.731 


0.777 


52 


105.2 


18.904 


13.045 


21 


45.0 


3.015 


0.908 


53 


105.8 


19.690 


13.813 


22 


46.0 


3.312 


1.050 


54 


106.4 


20.499 


14.613 


23 


47.0 


3.623 


1.205 


55 


107.0 


21.332 


15.446 


24 


48.0 


3.948 


1.372 


56 


106.4 


22.190 


16.314 


25 


49.0 


4.286 


1.551 


57 


105.8 


23.075 


1.217 


26 


50.6 


4.637 


1.744 


58 


105.2 


23.986 


18.155 


27 


52.2 


5.001 


1.950 


59 


104.6 


24.925 


19.131 


28 


53.8 


5.378 


2.170 


60 


104.0 


25.894 


20.145 


29 


55.4 


5.769 


2.404 


61 


98.0 


26.894 


21,199 


30 


57.0 


6.174 


2.652 


62 


80.0 


27.925 


22.293 


31 


58.8 


6.591 


2.916 


62.5 


20.0 


28.450 


22.850 



Table 2. Sample TBM Profile [Ref 5] 
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A profile consists of IR intensity as a function of time, nominal (maximum range 
trajectory) verti cal and horizontal ^ges from the launch point as functions of time, and 
maximum bum time, tmax (62.5 seconds in the example) The detected radiant intensity 
over time is compared to TBM profiles in a data base, and the best match is selected as the 
type of TBM being observed This selection process, called “typing”, is complex, and is 
also assumed to have been done correctly. A more computationally efficient (but equiva- 
lent) representation of the nominal downrange distance and altitude profiles can be found 
by fitting quartic polynomials to the data in the profile tables Use of the polynomial repre- 
sentation during calculations precludes table look-up and interpolation problems for non- 
integer times-of-flight. The profile distance polynomial has the form: 

dp =ao +a,t + a2t^ +ajt^ +a4t'* (1) 

and the profile height polynomial is: 

hp = b() + b,t + bjt^ + bjt^ + b^f* ( 2 ) 

where t is time of flight, and ai ,4 and bi .,4 are coefficients of the fourth-order polynomial 
These coefficients are not computed within this algorithm, but are assumed to be known, 
exact, correct, and available from the typing process; i.e., perfect a priori knowledge of 
the particular observed TBM’s nominal trajectory. Since real-life TBM’s do not fly the 
profile exactly, using the profile to determine the trajectory introduces an error into the 
algorithm The inherent error of inexact profile information is ignored in this analysis 
Assuming the profile to be exact does not invalidate the error analysis, however. 

A “loft” parameter, L, is used to account for trajectories above or below the nomi- 
nal profile. This very simple model for a ballistic missile trajectory adjusts the profile to 
give actual range. 

d = (l-1.5L)dp (3) 

and altitude: 

h = (l + L)hp (4) 

Loft varies over approximately ±0 25. Figure 6 is a plot of nominal (L = 0), lofted (L = 
+0.25) and depressed (L = -0.25) trajectories. 
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C. LINE OF SIGHT PROJECTION 

The least-squares iteration process begins with an initial estimate of the tactical 
parameters. In order to make an estimate for the initial TBM position (Xo, (|)o), the focal 
plane measurements must be changed into a LOS between the satellites and the TBM, 
which points to the position of the TBM on the globe The first step is to transform the 
satellite positions (gha, 6, R) into earth-centered rotating coordinates (XYZ), in which the 
X-axis extends through the prime meridian and the Z axis is aligned with earth’s spin axis, 
depicted in Figure 7 
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z 




Figure 7. Coordinate Frame Illustration 
The satellite position transformation from (gha, 6, R) to (XYZ) is; 



fx ^ 




^ cos(8 ^ ) cos(gha 


ysk 


= 


Rk cos(5^)sin(gha^) 






^ R^sin(5J j 



where, the subscript, s, denotes satellite position, and the sub-subscript refers to the k*^ 
observation. The focal plane coordinates, T|k and Pk, are first transformed to the satellite’s 
Up-East-North reference frame, (UEN), and then must also be transformed into (XYZ) 
The focal plane coordinates, can be visualized with Figure 8. 
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The transformation of polar (rit, 3k) coordinates into Cartesian (UEN) coordinates 



-C0S(T|J 

«k = -s»n(Pk)sin(Tik) 
-cos(3Jsin(T|J 



( 6 ) 



This is the unit line-of-sight (LOS) direction vector in (UEN) coordinates. Next transform 
into (XYZ) (refer to Figure 7): 





cos(ghaJ -sin(ghai) 


O’ 


cos(8k) 0 -sin(5J* 


Ck = 


sin(ghaj cos(ghaJ 


0 


0 1 0 




0 0 


1. 


sin(8^) 0 cos(8k) 



For the initial position estimate, it is sufficient to assume that the line of sight 
terminates on the surface of a spherical earth vidth an effective radius of reff = 6371 km. 
The Cartesian coordinates of the TBM are: 



= R, +p^e^ 



( 8 ) 



where the subscript, t, denotes TBM position, sub-subscript k refers to the k'*' observa- 
tion, and p is the length of the LOS vector as shown in Figure 7. The Law of Sines is one 
of several methods that may be used to calculate p. The geometry of the problem is 
shown in Figure 9. 




Figure 9. Initial Position Estimate Geometry 



16 



The angle opposite p is: 




( 9 ) 



With some trigonometric identities and algebra, the expression for p becomes: 



_ r,g.sin(xj 

sin(TiJ 

The longitude of this point is: 



r^g^ sin Ti^ + sin 




( 10 ) 



sin(Ti^ ) 




( 11 ) 



and the latitude is: 




( 12 ) 



D. INITUL ESTIMATES 

The core of the estimation process is a nonlinear least-squares iterative calculation 
of the tactical parameters. This method is based on a one-term (linear) Taylor expansion 
of the focal plane observations in terms of the six tactical parameters. This procedure 
requires an initial estimate of the tactical parameters. The accuracy of the initial estimate 
does not eifect the accuracy of the final estimate, but can effect the convergence time of 
the least-squares process. 

The initial estimate of the launch position and heading is based on the projected 
positions (<K and Xk) obtained fi'om the LOS observations. The lines-of-sight from each 
satellite usually do not intersect at a single point (even if they were simultaneous observa- 
tions), as shown in Figure 10. The initial latitude and longitude guesses are found by 
calculating the average position from the first LOS observation from each satellite. The 
assumption is made that if one satellite can see the TBM, then all can observe it. This 
artificiality is not true in the physical world. 
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Figure 10. LOS Projection Scatter 

In this example, three satellites observe the TBM, so the first LOS projected lati 
tudes and longitudes (from each satellite) are averaged: 

, (0) _ <l>| +<t>2 +<t>3 






A,| + X2 + ^3 



( 13 ) 



( 14 ) 



The superscript, denotes the initial estimate. The last LOS projected positions are also 
averaged to obtain the approximate final position observed: 

± _ <t>n-2 1 +<l>n 

Tlast “ ^ 



1 _ ^0-2 + I + 

^Ust ~ 2 



( 15 ) 

( 16 ) 
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Equations (13) to (16) assume that the first and last three observations are from three 
different satellites. If only two satellites observe the launch, only the first and last two 
observations would be averaged. 

The direction of the last position fi-om the first is the initial estimate for the head- 



= Y-tan,-'[(X,, -C^)cos(C>X<|,,, -C’] 



(17) 



The function tan 2 '' is the two-argument arctangent, used here because its range is ±tc. 
This function is used because a four-quadrant answer is required 

The first estimate for To is twenty seconds before the first observation: 

^ T/^=T, -20 (18) 

and the initial guess for L and ho are both zero: 

=0 (19) 

ho'°’ = 0 (20) 



E. EXPECTED POSITIOXS ALONG THE TRAJECTORY 

Given the six (estimated) tactical parameters at launch, the expected position of 
the TBM along its boost trajectory can be determined by applying the TBM profile to any 



observation time, and in particular, Tk, by first computing: 

tk=T, -T„ (21) 

Equations (1) through (4) then give the expected altitude (hk) and range (dk) fi'om the 
estimated launch point when t = tk 

dp, =ao +a,t^ +ajt^^ +ajt/ -^a^t/ (1) 

hp^ =b„+b,t, +b,t,' (2) 

d, =(l-1.5L)dp^ (3) 

h^=(l + L)hp^ (4) 
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Next, it is necessary to calculate the “earth-central angle” a quantity used in (23) and 
(24): 

©k = — (22) 

^eff 

0k is simply the angle between the launch point and the position of the TBM with the 
earth’s center as the vertex. 0 can be visualized with Figure 1 1 . 



V 



Earth ^ 





Launch 

position 



Figure 1 1 . Earth-central Angle 

Using 0k, spherical trigonometry gives the equations for the geodetic latitude; 



(|>^ = y-cos‘’[cos(0^)sin((|)o) + sin(0^)cos(<l)o)cos(ao)] 



and longitude: 



sin' 



,i| sin(0 Jsin(at,) 



(23) 



(24) 



cos(<l)J 

The altitude is simply the present height of the TBM added to the initial height at launch; 

alt, -ho+h^ (25) 

In addition to the geodetic coordinates, the Cartesian coordinates are also required Using 



re = 6378.137 km, and f = - 



1 



298.257 

the geocentric latitude of the TBM is calculated; 

= tan ^[(1 - f)^ tan(<j), )] 



the WGS-84 values for the flattening of the earth. 



(26) 
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as well as the local radius: 






( 27 ) 



■J(\ - f)’ cos^((t>\ ) sin^(<t)\ ) 

These values are then transformed into earth-centered, Cartesian coordinates of the TBM. 

Xt, = [rioci* cos((|)\ ) -t alt^ cos((t), )]cos(X^ ) (28) 

Yt, = ['"loci, cos(<j)\ ) + alt^ cos((|), )]sin(?., ) (29) 

sin(<l)\ ) + alt^ sin((l)^ ) (30) 



F. NONLINEAR LEAST SQUARES ESTIMATION 

Since the (polar) focal plane coordinates of the observations are nonlinear func- 
tions of the tactical parameters, a nonlinear least-squares process is required. This is 
achieved by a sequence of linear least-squares estimations. Each of these linear least 
squares estimations is based on a one-term (linear) Taylor series expansion of the obser- 
vations in terms of the tactical parameters A general requirement of ordinary linear least 
squares is that the noise contaminating the observations should be independent and have a 
common (preferably Gaussian) distribution. 

In this case, the “polar” focal plane observations ((3k, T]k) are transformed into 
“Cartesian-like” coordinates (Uk, Vk): 

u^ = -tan(Ti Jsin(P J (31) 

v^ = -tan(ri Jcos(P J (32) 

so that the coordinates have similar behavior with respect to errors and noise. Either r\ or 
sin(ri) could be used in place of tan(T|), but since t|< 10°, all three behave similarly. The 
chosen form corresponds to a gnomonic projection of the globe onto a plane tangent at 
the nadir point Figure 12 illustrates this 
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Figure 12. Focal Plane Transformation 

Next, it is necessary to determine what the focal plane coordinates would be for 
each observation if the TBM were actually located at the predicted coordinates ^ven in 
(28) to (30). 



The expected line-of-sight is computed: 





r, s, 


1 Ay, 


= 1 yr, - Ys, 


IAz,> 


Ur, 



(33) 



In order to determine the focal plane coordinates of azimuth, p, and elevation, t|. 



the LOS vector given in (33) must be rotated into the observing satellite’s local vertical 
coordinate frame, (UEN): 



E. 

InJ 



f cos(5^) 0 sin(8^)^|^ cos(gha^) sin(ghaj^) 
I 0 10 j -sin(gha^) cos(gha^) 

i^-sin(S^) 0 cos(6^)J^ 0 0 



AyJ 

AzJ 



(34) 



This is the transpose of the transformation matrices given in (7). 

Now u and v can be simply expressed in (UEN^ terms by noting that: 



Mh=-r== ( 35 ) 

cos(p,)^ (36) 
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( 37 ) 



(38) 



Expected focal plane coordinates, (u^ ), are then: 

E, 

u. = 

U. 

N, 

Expected focal plane coordinates are denoted with a “ ^ Taking the difference between 
actual observations and expected observations generates 5u and 6v (found on the left-hand 
sides of (41) and (42)): 

6u^ = - Uj (39) 

5v^ = (40) 

Now that the changes in the focal plane coordinates are known, they can be related to 
changes in the tactical parameters The one-term Taylor expansions of u and v in terms of 
the tactical parameters are; 



5u = — 5T„ 
dv 

Sv = ^6T, 

ar. 



5L- 



5L 



au 

«t>o 

dv 



5(t)o- 



5«i)o 






5:^0 



au 

av 



5h„ 



8h„ 



au 

h — 

da^ 

dw 






8a„ 



(41) 



(42) 



aL a4>o " dk, ° ah„ " aa„ 

These equations, written in vector-matrix format, are the linear least-squares 
model that forms the basis of the tactical parameter estimation process: 



^ 5ui 
kSv^ 



f 


5Uj. 






du. 


1 

du^ Y 


1-1 


dL 




dx. 


dh. 


da^ 1 






^k 


cVk 




^k i 


UtTo 




^0 


dX, 


dho 


daJl 



5L 

5(t>o 

5?.o 

5ho 



(43) 



There are two rows of the center matrix in (43) for each of the k = (1 to n) observations 
for a total of 2n rows. This matrix of partials is denoted the “A” matrix, and represents 
changes in the focal plane coordinates with respect to changes in tactical parameters. 
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The differences in the focal plane coordinates, Su and 5v, are used to generate 
adjustments to the tactical parameters by solving (43) for the tactical parameter variations; 



6L 

&|)o 

She 



■K«i '*’£:) 



(44) 



The middle term in (44) is the pseudoinverse of A, and is used because A is an (2n x 6) 
matrix, and is generally not square. The changes to the tactical parameters, the left-hand 
matrix in (44), are added to the initial estimate, and the process is repeated until all the 
components of the left-hand vector in (44) are sufficiently small, or the convergence crite- 
rion is met. The result is the correct tactical parameters for the observed TBM launch. 

The matrix of partial derivatives. A, is the key to the whole iteration process. The 
development of this matrix involves multiple use of the Chain Rule for derivatives. The 
partial s in (41) can be expanded to yield: 



5u 


^ ^ 5u du dh 

~ d^dl^ ' dk dT, ^ dh dY, 


(45) 


5u 


^ ^ du dk du dh 


(46) 




^ 5L dk dL dh dL 


c*u 


^ ^ da dk 0u ?h 


(47) 


^0 


0<1> Scj)^ dk 0h a|>(, 


du 


0u ^ da dk cri ^h 


(48) 




c5(j) dk^ dk dk^ 5h dk^ 


du 


0u 5(|) da dk 5u 5h 


(49) 


d\ 




du 


c^o(|)^^a^5u6h 


(50) 


da^ 


5(1) 5a„ dk 5ap 5h da^ 
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In a similar manner, the partials of (42) are expanded: 

dv dv dv dk dv 

^ ^ ^ a 


(51) 


dv dv dv dk dv dh 


(52) 


dv dv d^ ^dv dk ^dv dh 

r4o ^ ^0 ^0 ^ ^0 


(53) 


dv f\|) dv dk dv ^ 

^ a ^ 


(54) 


dv ^ t\j) dv dk dv dh 

flhp t\t> ^0 ^ ^^0 


(55) 


dv dv d^ ^dv dk ^dv 

da.Q oj) da^, dk da^ rlh da^ 


(56) 


Patterns and structure in these twelve equations are more easily recognized when 



these equations are written in matrix form: 



^au^ 


r f\t) 


dk 


ah ■ 


! axo 


1 axo 


ax„ 


axo 


au 


at) 


dk 


ah 




dL 


dL 


dL ( 


au 




dk 


dh 


ax{)o 


(^0 


a^o 


at>o 


au 


at) 


dk 


ah 


1^ 


dk. 


dk. 


dk. 


au 


a|) 


dk 


1 


~^0 


aho 


dh. 


aip 


du 


axt) 


dk 


ah 


vaao> 


da. 


da. 


axo 
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dk 


dh ' 




1 ^0 


ar„ 


dT, 


5v 


1 


dX 


dh 




dL 


dL 


dL ( 


dv 




dX 


da 




\ 


at)o 




dv 


“ d^ 


dX 


ah 




\ ^^0 


dX, 




dv 


d^ 


dX 


ah 1 






dh. 


ahp 


dv 


ad) 


dX 


ah 


\daj 


da. 


da. 


da. 



Note here that (57) and (58) have identical structure. The left-hand sides are the deriva- 
tives of the focal plane coordinates v^th respect to the tactical parameters, the vectors on 
the right are the derivatives of the corresponding focal plane coordinates with respect to 
TBM position (latitude, longitude, and height), and the center matrix describes the deriva- 
tives of the position with respect the the tactical parameters. This effectively separates the 
change in the focal plane coordinates with respect to tactical parameters into two parts: a 
matrix that describes changes in position due to trajectory, and a vector that describes the 
relation between focal plane observations and TBM position. 

The matrix is given a name, Ai: 



■ ^ 


dX 


ah ■ 


dY, 


aTo 


dT, 


d^ 


dX 


dh 


a 


dL 


dL 


atj) 


dx 


ah 




(^0 


Wo 


atj) 


dX 


ah 




dx. 


w 


a(|) 


dX 


dh 




dh. 


dh. 


d^ 


dx 


ah 


da. 


da. 


da. 
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so that (57) and (58) can be rewritten: 











1 da 




' aL 




^au^ 








1 ^0 




au 


j da 


= Aj 








ail 






xah> 


‘ aho 




‘ au 








r av ^ 




: 




av 








^ d\'' 


C\’ 








= A, 


dv 


dv 




CK 






C\ 


cv 










d\' 









( 60 ) 



(61) 



The elements of Ai can be found qualitatively This is done in Appendix A. and 
the results are; 



od cos(ao) 

^0 fcff 

1.5dpcos(ao) 



1 

0 

0 

dsin(a„) 



dd sin(gp) 
cT^ cos((|)o, 

1.5dp sin(g(,) 
r^cos((|)„) 

0 0 

1 0 

0 1 

dcos(ao) Q 

r^cos((|)o) 



(62) 
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The values of the elements in (62) are approximations, but approximate magni- 
tudes and correct signs are all that are required for the iteration process to converge. The 
benefit of making this simplification for Ai is the reduction in computation time The 
complete and exact formulae can be extrapolated fi-om Appendix A 

Similarly, vectors on the right-hand sides of (57) and (58) can be separated into 
simpler components: 



^ _ 


0U 


dx 


du 




au dz 




d\ 


dy 


at) a at) 


du 


du 


dx 


du 


dy 


du dz 


a ' 


dx 


a' 


^ dy 


dX 


^ a a 


<5u 




dx du 


dy ^ 


a dz 


0h 


'a 


ah dy 


ch 


dz ah 


d\ 


dv 


dx 


d\ 


dy 


dw dz 


at> 


dx 


at) dy 


at) az at) 


dv 


_ ^ 


dx 


C\' 


ay 


d\ dz 


Ijk' 


dx 


a' 




dX 


dz dX 


dv 


_ dv 


dx 


dv 


dy 


c\ dz 


a ' 


dx 


a' 


^y 




dz ah 



As before, these equations can be written in matrix form: 







dx 




a 








at> 


at) 


at) 


6X 


a 




dx 


dy 


a 


a 


a 




a 


dX 


a 


dy ' 


a 


1 

1 


dx 


h. 


az 


— \ 


laj 


i 


a 


a 


a 


'^dzj 


^ dv^ 


1 


dx 


dy 


a 


fdv^. 






at) 


at) 


di 


dx ' 


dv 




dx 


dy 


dz 


av , 


a 




a 


a 


a 


dy 


dv 




dx 


dy 


dz 




la. 




a 


a 


a 


— 

^az>' 



(63) 

(64) 

(65) 

( 66 ) 

(67) 

( 68 ) 



(69) 



( 70 ) 
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The center matrix is renamed A2: 



’ax 




dz' 


a«j> 






ax 


ay 


dz 


a 


dx 


dx 


ax 


ay 


dz 


ah 


ah 


ah 






and its elements can be approximated by: 

’-(reff + h)sin(4>)cos(X) -(r^ + h)sin(4>)sin(X) (r^ + h)cos(4>)‘ 
Aj = -(r^ + h)cos(<t>)sin(X) (r^ + h)cos(4>)cos(A.) 0 

cos((|))cos(X,) cos((t>)sin(X) sin(4>) 

The vectors on the right-hand side of (69) and (70) can be further broken down: 

(?U 

VcN. 







au 


dE 


aN' 


ax 




dx 


dx 


dx 


au 




au 


dE 


aN 


ay 




ay 


dy 


ay 


CU 






dE 


aN 


\dz^ 




. dz 


dz 


dz . 


^ d\^ 




'au 




aN" 


chi 




dx 


dx 


dx 


av 




au 


dE 


aN 








dy 


dy 


dv 




au 


dE 


aN 


vaz> 




. dz 


dz 


dz . 



The center matrices in (73) and (74) are identical, and renamed A3: 

aN' 



(71) 



(72) 



(73) 



(74) 



A^ = 



^ ^ 

ax ax ax 
^ ^ ^ 
ay ay ay 
au ^ ^ 

,dz dz dz 

A3, th transformation from (XYZ) to (UEN) has been derived previously in (7): 



A 3 = 



cos(gha^) -sin(gha^) 0 
sin(ghaj cos(ghaJ 0 
0 0 1 



cos(5^) 0 -sin(5^) 
0 1 0 
sin(5J 0 cos(5^) 



(75) 



(76) 



29 



The final step is to determine the elements in the vectors on the right-hand sides of 
(73) and (74). These are easily computed from (37) and (38): 



^ = A 
au U' 

au _ j_ 
ai ~~u 

A. 

aN ' 
av _ N 

^ = 0 
aE 



- = o 



dv 

lis' 

Now (57) and (58) can be rewritten: 

^ du ' 

au 

dL 






aho 

au 



u' 



= A,AjA3 



Kcxxo^ 



(77) 

(78) 

(79) 

(80) 
(81) 
(82) 



(83) 
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f dv'' 

I ^ 

dv 

'K 

dv 

I dv 

1^ 

I dv 
I aio 



I 

0 , 
-—I 

uj 



\da^^ 



(84) 



Equation (83) is used to compute the odd rows (k = odd) of A in (43), and (84) is 
used to compute the even rows (k = even) of A in (43). 

Focal plane coordinates, (3k, tlk), are transformed into Cartesian-like focal plane 
coordinates, (Uk, Vk) Cl: (Circled numbers refer to steps illustrated in Figure 13 ) Satellite 
positions and initial estimates of tactical parameters are used to generate expected 
(theoretical) focal plane coordinates, (u^ , v^^ ) S®, and the A matrix (which is denoted 

by in Figure 13) and its pseudoinverse S . The differ ences of the 

rXTo,L,4>o,Xo,ho,ao) 



focal plane coordinates, (5uk, 5vk), (observed minus theoretical) 'S' are used to generate 
adjustments to the tactical parameter estimates by computing 





^5u,^ 




6v, 


SL ; 


8u2 


5<l>n -r 1 tI 


Svj 


i ;=[AU]'aX 




i 5h„ ^ 






5u„ 


1 





If all the elements of the left-hand vector (adjustments to the tactical parameters) are not 
sufficiently small, the current values 0- 1 ) of the tactical parameters are updated ® to the 

f: 
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(86) 



1 ° 1 


(T ^ 




L 


L 


: 8L 


<l>o; 


4>0 










^0 


ho 


i»l>„ 




UoJ 





If the iteration process has converged, the iteration process is terminated and the current 
values of the tactical parameters are the best estimates that this process can generate. This 
entire process may be visualized by a flowchart diagram drawn in Figure 13. 




Figure 13 Flowchart Diagram of Iteration Process 



G. BURNOUTTIMEESTLMATION 



The estimation of burnout time is based on the TBM profile maximum bum time 
(tmaT), the last observation time (Tiaa), and the next potential observation time (Tort), had it 
occurred The maximum bum time according to the profile is; 

T_=To+t_ (87) 

Two cases can occur: 

(a) ifT:M.x>T„e^ then 

T>.=T^, (88) 
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(b) if T,„„ < Tncxi then 

T^=T^'|(T„-T^> (89) 

A figure and example are given below as a demonstration: 




Figure 14 Burnout Time Estimation 

For the sample data given in Table 1, using the profile given in Table 2, t^ax = 62 5 
seconds, Tia^ = 155.44, and the next potential observation is T„cxi = 159.36. With To = 
109.36, the maximum bum time is = 171.86, so case (a) applies and burnout time is 
estimated at Tbo = 157.40. 

H. STATE VECTOR GENERATION 

The state vector completely defines the TBM’s position and velocity, and has six 
elements With the tactical parameters and burnout time estimates, generating the state 
vector at burnout is done by evaluating the position and velocity equations at the burnout 
time: 

t.o=T,„-T, (90) 
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Substitute tbo for t in (1) and (2), and use (3), (4), and (22) through (25) to get 4)bo, Xbo, 
and altboi 



+a,t^ +a,tfc„^ 


+ a^t J 


(1) 


hp, =b„+b,t^ 


. +bjt J +b3t J 


+ b,tj 


(2) 


dbo=(l 


-1 5L)dp^ 




( 3 ) 


hbo=(l + L)hp, 




( 4 ) 








(22) 


I 


eff 






^ - cos ' [cos(0^ )sin((j> 0 ) + sin(9 ^ )cos(4>o )co%{a ^ )] 


(23) 


^bo = Xo ^sin'' 


sin(0bo)sin(ao) 

COS(<|)^„) 




(24) 


altbo ^ 


“ b(j h^^ 




(25) 



The burnout velocity is expressed in terms of speed (Vbo), flight path angle (ybo), 
and heading (Obo): 



Vb.. = Vd- + h- 

Ybc =»aii- 



'll 



(91) 

(92) 



^ ^ / cos(itio)sin(ao)' ; 

I cos(0b„) J 

where d and b are the time derivatives of d and h 

d=« 

a 



h = 






(93) 



(94) 

(95) 
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The state vector at burnout, Tbo, is 



fK] 

^bo 



Ybc 

L CALCULATION OF IMPACT POSITION AND TIME 

{ The equations and methods in this section are adapted directly from Chapter 6 
of Bate, et ai) [Ref 6] The ballistic trajectory is modeled in three phases: powered 
flight, which is the portion that DSP observes, free-flight, which is a portion of an elliptical 
orbit, and re-entry, which is the portion from where atmospheric drag becomes significant 
until missile impact. The powered-flight phase is modeled v^nth the TBM profile polyno- 
mials presented earlier. The ellipse traced during ffee-flight is simulated using inertial 
two-body mechanics Atmospheric drag effects during the re-entry phase are not specifi- 
cally calculated in this model, but are somewhat accounted for by assuming that the dis- 
tance traveled over the earth from re-entry to impact is the same as from launch to 
burnout This is the same as assuming that the earth-central angles are equal 

6re=0bo (97) 

The speed of the TBM during re-entry is assumed to be unaffected, however These 
approximations are recognized artificialities, but are simpler than calculating ballistic 
coefficients of various TBM’s and modeling atmospheric density, and more accurate than 
assuming the missile remains on its elliptical path until impact. Figure 15 illustrates the 
geometry and some of the quantities used in the equations. 
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Figure 15. Ballistic Trajectory [Ref 6] 

The equations contained in Chapter 6 of Ref 6 concern ballistic trajectories and 
are based upon inertial quantities and a spherical earth model, so the geocentric position 
vector and inertial speed at burnout are required. Geocentric latitude is generated from 
(26): 

4>\„ = tan ’[(l-f)^tan(<|),„)] (26) 

The earth’s radius at burnout is calculated with (27): 

, - p,) 

'locii^ / ; ; ; 

V(1 - f)^ cos^(f bo ) + sin^ (4>'b„ ) 

These values are then transformed into earth-centered, Cartesian coordinates of the TBM: 
cos(4>'fc„ ) + alt^ cos(<|)^)]cos(>.^) (28) 

Vt, =[rk,«^„cos((t>*^) + aIt^cos(<|)^)]sin(?L^) 

= W ) + alt^ sin((j>^) 
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(29) 

(30) 



which define the geocentric TBM position at burnout; 

I ( 8 ) 

In an inertial reference frame, the TBM has an initial eastward velocity at launch 
due to the earth’s rotation, Vo. This velocity can be expressed as: 

Vo = r,^,^ti>,cos(4>;j (98) 



where tOe = 15°/hour, the rotation rate of the earth The burnout velocity can be broken 
down into three inertial (UEN) components. In the inertial Irame, the initial eastward 
velocity is added to the eastward component 

V. =V^sin(Y^) (99) 

"Ve = V^cos(y^)sin(a^) + v„ (100) 

= V^cos(7b„)cos(a^) (101) 

Inertial quantities are denoted with a bar, “ ~ The magnitude of this vector is the iner- 
tial speed, 



V,. = Vv," - + V,’ 

The inertial flight path angle and heading may now be calculated: 




( 102 ) 



(103) 



= tan i ! 



(104) 



Now that the inertial quantities are known, a non-dimensional parameter, Q, is defined as 
the squared ratio of the inertial speed of the TBM to circular orbit speed at that position: 



Q = 



It 



(105) 



where p = 398601.2 — the gravitational parameter for the earth, and rbo = irbo|, the 
sec 

magnitude of the TBM position vector. 
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The eccentricity of the ballistic “orbit” is: 



e = v/^+Q(Q-2)cos"(7bo) 



(106) 



The free-flight “earth central angle”, T, is defined by Bate, et al., as that earth- 
central angle that the missile traverses between burnout and re-entry (see Figure 1 5). This 
definition is modified by adding the angle traversed during reentry, 0„ to the BMW de- 
fined angle, 4^ This implements the assumption that the TBM’s tr^ectory is modified by 
drag. The new definition for the angle, T, is: 



T = 2cos 



, p-Qcos^(Y^o)j ^ 



(107) 



The eccentric anomaly, E, and the semi-major axis, a, of the ballistic trajectory 







(108) 



a 



2-Q 



(109) 



Assuming the speed of the TBM is unaffected by drag during re-entry, the time of 
free flight (burnout to impact) is defined as: 



,=2 [7r-E + esin(E)] 



( 110 ) 



the time of impact fi-om launch is 
and the time of day of impact is: 



tun tbo tfr 

Tim ~ To + tim 



Using the Law of Cosines from spherical trigonometry, latitude at impact, is: 
4>L = sin''[sin(4i;,,)cos(4^)^cos(27i-a^)] 

This can be transformed back into the WGS-84 latitude: 




( 111 ) 

( 112 ) 

(113) 



(114) 
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This method is not an exact transformation, but it is a close approximation that does not 
affect the error analysis. A better model should be used if more precise latitude of impact 
is desired 

Using the Law of Cosines again, and after some algebra, the longitude traversed, 

AX, is: 



Better models for the ballistic trajectory could possibly be employed. The goal of 
this analysis is not to precisely determine the actual impact zone, but to determine the 
effects of error sources upon the result Thus, the effects found using this model should 
be applicable to other ballistic trajectory models. 




( 115 ) 



So, the impact longitude, Xan, is: 



^nn ~ ^bo 



( 116 ) 
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IV. ERROR SOURCES 



Noise (error) is present in every aspect of the algorithm presented in Chapter 3 
Every measurement, constant, coefficient, model, equation, approximation, and 
assumption impart some finite quantity of uncertainty upon the result. Even if a quantity is 
exact to its limit of accuracy, the simple fact that it is represented by a finite number of 
digjts limits the precision of that quantity For example, if a ruler has one-sbcteenth inch 
gradations, the accuracy of any measurements made with the ruler is ±one-thirty-second of 
an inch This means that there are no “perfect” values in the algorithm, and each value 
used is a source of uncertainty, termed an error source 

Each error source usually has one or more underlying causes A complete error 
analysis would break each error source down to its fundamental level, and model each 
level correctly An analogy for this is the layers in an onion: The onion can be viewed 
externally as a whole, but can also be peeled, layer by layer, to reveal deeper, underlying 
matter that support the external skin In this section, some of the underlying causes of the 
overall error sources are identified, but in these analyses, only the overall error magnitudes 
will be modeled and analyzed. 

The obvious origin of any errors present is the observational data since it is the 
starting point for the algorithm These data are physical measurements of three types 
time, focal plane measurements, and satellite position and orientation 

Time errors can be caused simply by having more than one clock being referenced 
as a source of time measurement Imperfect synchronization between clocks’ time and 
time passage rate are obvious error sources Time delays caused by radio transmission of 
data due to distance, atmospheric refi’action, and relative motion Doppler effects may add 
another time error. In the extreme, the source of time error can be traced all the way back 
to our measurement of time passage with respect to the celestial sphere, which is not 
fixed. The magnitude of time errors is relatively small, and getting smaller as time 
measurement improvements are being made continually 
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The largest source of time error enters the algorithm at the burnout time estimation 
phase The 10-second sampling rate limits the accuracy of the estimate to :t 5 seconds for 
single satellite observations. The accuracy of the profile’s value for tma also comes into 
play in the estimate The burnout time is crucial to determining the ballistic trajectory, 
since the TBM is undergoing its greatest acceleration changes during the final seconds of 
powered flight. Small errors in the estimate of burnout time thus result in larger errors in 
the estimate of impact time and position. 

LOS measurement errors can arise fi’om many noise sources, mainly of two types: 
attitude uncertainties and IR radiation measurement errors. The attitude of the spacecraft 
needs to be precisely known (and it simplifies the process if it is very stable). A 28 
microradian error in pitch or roll pointing accuracy of a geosynchronous satellite equates 
to a 1 kilometer error on the surface of the earth, at least 35786 kilometers away. Any 
attitude control system inaccuracy effects are amplified by the geosynchronous altitude. 
Vibrational noise fi-om spinning momentum wheels, uncertain knowledge of the mass 
properties of the yaw-spinning satellite, star catalog and star sensor measurement errors, 
thruster misalignments and disturbance torques, tachometer errors, and mechanical 
misalignments aD contribute to the total attitude uncertainty The knowledge of the 
telescope boresight axis alignment with the satellite’s reference fi-ame is defined by the 
design and manufacturing of the DSP satellite, and changes slightly with thermal 
variations. 

The ER radiation measurements of intensity and angle of arrival (AOA) also have 
several underlying error sources Locations of the individual photoelectric cells on the 
focal plane array are recorded in what is termed the ‘Tocal Plane Vector Table” (FPVT). 
The positions of the PEC’s change as the satellite heats and cools with varying sun- 
satellite orientations, causing a warping of the focal plane, but the FPVT does not account 
for the changes real-time. In addition, detector noise, refi-action, and IR attenuation due 
to clouds and water vapor all add to the total measurement error. 

Satellite position measurements fi-om AFSCN are used to update orbit ephemeris 
data. The measurements are, of course, inexact, but those errors are compounded over 
time because the updates occur only once weekly The ephemeris data is propagated to 
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determine satellite position during the week between updates DSP is a geosynchronous 
satellite, so hs orbit does not change greatly in that period, but errors in the propagation 
model can cause the predicted position to be different from reality. 

Keeping all of these underlying components in mind, the errors are modeled in the 
tactical parameter estimation algorithm The real-world error component magnitudes, bias 
and random distributions may be different from those simulated, but the overall effects 
manifest themselves in a manner close to the error models. It should be possible to 
extrapolate the results obtained in this study to different errors by properly scaling the 
different magnitudes and distributions. It must be emphasized that since the goal of this 
study is not to exactly determine the tactical parameters at launch, state vector at burnout, 
or impact time and position The purpose is to determine the contribution of the error 
sources upon these values, and this can be done with approximate models for the errors 
If the intended goal is to calibrate the system to eliminate bias errors, then the sources of 
the errors must be determined more exactly to determine what is observable (measurable) 
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V. NUMERICAL ANALYSIS 



The algorithm described in Chapter III was programmed into MATLAB™ Three 
error sources (time, satellite position, and LOS) were simulated and added to the 
observational data, first separately, and then combined. The effects of the errors upon the 
results are analyzed at three points: launch position, burnout position, and impact 

position. Five Middle Eastern capital cities were chosen as fictitious launch sites to 
determine the effects of TBM launch position upon the accuracy of the results 

The model for time error is a random uniform distribution between 0 and 1 
milliseconds. For satellite position, a random normal distribution with zero mean and 
standard deviation of 200 meters was added to each of the components (R, gha, 5). This 
is approximately a one standard deviation sphere with radius 346 meters. The LOS error 
was modeled as a random normal distribution with a standard deviation of 5 microradians, 
added to both (5 and rj components Histograms of sample error distributions are shown 
in the following three figures. 




0 02 04 0.6 08 1 



Tinrie error (seconds) ^ ^q-3 

Figure 16. Histogram of Time Error - Uniform Distribution 
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Satellite Position Error (meters) 

Figure 17, Histogram of Satellite Position Error - Normal Distribution 



0.7 r 
0.6 - 




Focal Plane Error (radians) ^ ^g-s 

Figure 18. Histogram of LOS Error - Normal Distribution 
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Errors added to the observations had the effect of creating an ellipse of incorrect 
data points around the true position, termed an “error ellipse” The l-o ellipses shown in 
the next figures are determined using statistical methods. [Ref 7] For the launch ellipse, a 
matrix of two columns was collected during the simulation: 






lpts = 




(117) 






where Ipts = launch points, the matrix of launch longitudes, Xo, and latitude, (j)o. The first 
row is the generated true launch position, the second row is the calculated launch position 
with no error added to the observational data, and the remaining 999 rows are the 
calculated launch positions with errors added to the observational data. The first two 
rows are error-free, and are removed to produce a matrix of only “corrupt” positions To 
determine the ellipse dimensions, a and b (semi-major and semi-minor axes lengths), the 
eigenvectors and values of the covariance of this matrix were calculated. Twice the 
square root of the eigenvalue diagonal elements produces a and b: 

2 .peigenvalue(l.l) 0 1 

I 0 ^eigenvalue(2,2) 




and the eigenvector is the direction cosine matrix for rotating the ellipses from the primary 
axes. The impact ellipses are calculated similarly, with the initial data matrix being 
composed of impact longitudes and latitudes, instead of launch The ellipses can then be 
plotted using the ellipse equation; 




(119) 



and letting x vary between ±a, and letting y be the dependent variable. Multiplying the 
resultant (x, y) coordinates by the eigenvector matrix rotates them to the correct 
orientations. The semi-axes are in units of radians, since the positional values are in 
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radians. Thus, the (x, y) plots produced are the longitude and latitude of the points, once 
the values are converted to degrees. 

For the burnout position, the same method is used, but the data point matrix has 
three columns, the third being the altitude of the TBM at burnout. The covariance method 
works as well in three dimensions, producing a, b, and c, the three semi-axes for the 
ellipsoid formed. The coordinates are calculated in the same manner, using the equation 
for an ellipsoid; 






( 120 ) 



and, again, rotating by pre-multiplying by the eigenvector matrix. 

The following figures are representative of all the simulated cases, but this 
particular case is a 300-km TBM launched from Baghdad heading 0°, with the combined 
erorrors 



Launch Position 




Figure 19. Launch Position Error Ellipse 

The ellipse drawn on top of the data points encompasses the data points within one 
standard deviation of the mean, assuming the distribution is normal. Since the time error 
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has a uniform distribution, 1-a is not the expected 66% for time error cases. This method 
is used anyway for lack of a better one and to maintain a common baseline for 
comparison. As it turns out, the time error contribution is very small, and thus does not 
effect the combined error case distribution. 

The position at burnout is plotted in three dimensions, using altitude as the third 
coordinate to form an ellipsoidal error volume. 

Position at Burnout 




Figure 20. Position at Burnout 

The three ellipses show the outline of the ellipsoidal volume. It is difficult to 
perceive the ellipse orientations, but they are mutually orthogonal, (necessarily, since they 
correspond to the eigenvectors of distinct eigenvalues). 
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The impact ellipse is similar but larger than the launch position ellipse, due to the 
propagation of velocity errors from burnout to impact and the error in the burnout time 
Impact Position 




Figure 2 1 , Impact Position Error Ellipse 

The three previous figures illustrate the error ellipses(oids), but showing all the 
data obtained in this format is unwieldy. The sizes (areas and volumes) of the positional 
errors are the important quantities for identifying which error source has the largest effect. 
Thus, the data have been processed to determine these quantities, and are plotted next. 

The data for the 300-km TBM was collated by city, heading, and error source. 
The ellipse(oid) dimensions calculated must be equated to the actual distance on the 
ground to determine the area of the ellipses. To do this, the angles must be in radians, and 
are multiplied by riocai The distance between lines of longitude shrinks as the latitude 
moves away from the equator, so the longitudinal distance must be multiplied by the 
cosine of the geocentric latitude, (j)' Thus the equation for ellipse area becomes: 

area = ;rr,„^,,^abcos((|)') (121) 
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The ellipsoid volume is not as simple, so a parallelepiped with dimensions, a x b x 
c, is used to approximate it; 

volume = abc (122) 

These areas and volumes are plotted in the following figures, separated by error 
source and heading. Each launch point is represented by a different color: 

Aden = blue 
Baghdad = black 
Damascus = green 
Riyadh = red 
Tehran = cyan 
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Burnout Ellipsoid Volume (knV'3) 
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Figure 27. Satellite Position Error Effects Upon Impact Ellipse Area 
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Burnout Ellipsoid Voluine (km^S) 
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Burnout Ellipsoid Volurrie (km^3) Launch Ellipse Area (km*2) 
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Figure 3 1 Combined Error Effects Upon Launch Ellipse Area 
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Heading (degrees) 

Figure 33. Combined Error Effects Upon Impact Error Ellipse 
It is immediately apparent from these plots that the launch site has little effect upon 
the size of the error ellipses and ellipsoids It is also obvious that the size is dependent 
upon TBM heading and observation geometry. The relative magnitudes of the error 
effects are not as obvious, so a plot of four error ellipses is shown in the next figure, with 
each ellipse formed by a different error source. 
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Launch Position 




Longitude (degrees) 



Figure 34. Launch Error Ellipses for Each Error Source 
From this plot it is easily recognized that the errors can be ranked in order of 

effect: 

1 . focal plane errors 
2 satellite position errors 
3. time 

All the error sources together produce the largest ellipse, as would be expected from the 
superposition principle. 



This is further illustrated by Table 3. Each error source is a row; “time” is the 
time error runs, “satpos” is the satellite position error runs, “LOS” is line of sight, and 
“combin.” is the combined errors runs 



ERROR 


LAUNCH ELLIPSE AREA 
(km^) 


BURNOUT ELLIPSOID 
VOLUME (km^) 


IMPACT ELLIPSE AREA 
(km“) 


min 


mean 


max 


min 


mean 


max 


min 


mean 


max 


time 


3.26e-8 


2.59e-7 


1.22e-6 


2.9e-ll 


3.6e-I0 


1.33e-9 


I.57e-6 


2.6Ie-5 


1.41e-4 


satpos 


3.77e-3 


5.27e-3 


8.68e-3 


6.04e-4 


I.I4e-3 


2.32e-3 


6.l2e-3 


7.92e-2 


3.27e-I 


LOS 


2.97e-2 


8.45e-2 


1.83e-l 


2.37e-2 


9.33e-2 


2.09e-I 


3.05e0 


I.02el 


I.95el 


combin. 


3.91e-2 


l.Ole-1 


2.97e-l 


4.44e-2 


I.25e-1 


4.50e-l 


3.35e0 


1.38el 


4.06el 



Table 3 Error Ellipse Areas and Ellipsoid Volumes 
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Ellipsoid Volume 

Ellipse Area (km^) (km») EUipse Area (km ) 



This is also shown graphically in the following plots: 




Figure 35 Launch Ellipse Area Comparison 




Figure 36. Burnout Ellipsoid Volume Comparison 




Figure 37. Impact Ellipse Area Comparison 
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The sum of the error sources’ positional error sizes is less than the combined error 
source’s positional error size. This implies a nonlinear interaction of error effects with 
each other. This is a reasonable result considering the complexity of the simulated system, 
and noting that most real systems exhibit nonlinear behavior. 

Taking advantage of the ease of modifying MATLAB™ code, various changes to 
the present DSP system model were put into simulation to determine the effects upon the 
accuracy of the results. For every case, Baghdad is the launch site and the observational 
data were modified by combined error sources. 

The first case is the control case with nominal system parameters. It is listed as 
“Bagcom” to represent “Baghdad combined errors”, and is used as a baseline case for 
comparison. The second case, “Synchr.” shows the effects of synchronizing the spins of 
the satellites so that the satellites scan a single area of interest within one second of each 
other For the third case, “Molniya”, the satellite at 70° GHA was elevated to 63 .4° decli- 
nation, to simulate a Molniya + geostationary viewing geometry. The remaining cases 
were simulating faster scan rates from 10 seconds down to 2.5 seconds. The numerical 
results are shown in Table 4 



ERROR 


LAUNCH ELLIPSE AREA 
(km^) 


BURNOUT ELLIPSOID 
VOLUME Ocm^) 


IMPACT ELLIPSE AREA 
(Icm^) 


min 


mean 


max 


min 


max 


max 


min 


mean 


max 


Bagcom 


4.216-2 


l.Ole-1 


2.746-1 


4.46e-2 


1.096-1 


2.60e-l 


3.83eO 


1.40el 


4.04el 


Svnchr. 


4.406-2 


7.076-2 


l.lOe-1 


2.096-2 


4.896-2 


8.616-2 


3.94e0 


7.82e0 


1.52eO 


Molniya 


4.396-2 


7.586-2 


1.256-1 


2.616-2 


5.576-2 


1.046-1 


2.91eO 


7.96e0 


1.85el 


10 s sc 


4.216-2 


l.Ole-1 


2.746-1 


4.466-2 


1.096-1 


2.60e-l 


3.83eO 


1.40el 


4.04el 


7.5 s sc 


2.626-2 


7.446-2 


1.956-1 


2.53e-2 


4.55e-2 


1.096-1 


4.03e0 


8.72e0 


2.08el 


5 s sc 


2.486-2 


5.026-2 


1.526-1 


1.236-2 


2.526-2 


6.52e-2 


2.52e0 


5.84e0 


1.64el 


2.5 s sc 


1.186-2 


2.36e-2 


3.486-2 


5.056-3 


9.216-3 


1.316-2 


1.03e0 


2.95e0 


4.58e0 



Table 4. Various Case Comparison 



The effects of the modifications are more easily observed graphically . 
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Ellisoid Volatile (km’) 




Mean Launch Ellipse Areas 



Figure 38 Comparison of Mean Launch Ellipse Areas 




Mean Burnout Ellipsoid Volume 



Figure 39. Comparison of Mean Burnout Ellipsoid Volumes 
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□ Baghdad Combined 
Gi Synchronized Spins 

□ Moliiiya 



Mean Impact Ellipse Area 



Figure 40 Comparison of Mean Impact Ellipse Areas 

Both modifications have the effect of decreasing the areas and volume by about a 
third This is the expected result for the Molniya case, since the viewing geometry is more 
three-dimensional with two satellites on the equator and the middle satellite near Molniyan 
apogee. 

The synchronized spin results were a surprise, however. The effect upon the 
launch ellipse area is fine, but since the burnout time estimation should have been less 
accurate, the burnout ellipsoid volume and impact ellipse area were expected to remain 
unchanged at best. Obtaining simultaneous observations would have obvious advantages, 
since they could be processed differently and be used to determine the position exactly for 
discrete instants in time during the boost trajectory. However, these observations were 
not processed differently, and were not simultaneous — only close together (within one 
second of each other). Still, the result may be a manifestation of some mathematical 
process that causes the increase in accuracy. 

The following plots show the effects of increasing the scan rate, or alternately, 
decreasing the time between observations This allows more observations to be made 
during the boost phase. 
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Ellipsoid Volume (km’) 




□ 10 s scan 
0 7.5 s scan 

□ 5 s scan 

□ 2.5 s scan 



Mean Launch Ellipse Areas 



Figure 41 . Scan Rate Effects on Launch Ellipse Area 




□ 10 s scan 
H 7.5 s scan 

□ 5.0 s scan 

□ 2.5 s scan 



Mean Burnout Ellipsoid Volumes 



Figure 42, Scan Rate Effects on Burnout Ellipsoid Volume 
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□ 10 s scan 
O 7.5 s scan 

□ 5.0 s scan 

□ 2.5 s scan 



Mean Impact Ellipse Areas 



Figure 43 . Scan Rate Effects on Impact Ellipse Area 



The effects seem to be linear decreases upon the launch and impact ellipse areas, 
and quadratic decreases on the burnout ellipsoids. This result is in perfect agreement with 
intuitive predictions. 

In summary, the numerical results are generally in accordance with expected re- 
sults. As a qualitative check of the quantitative results, the next chapter is devoted to a 
qualitative analysis of a simplified case for comparison. 
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VI. QUALITATIVE ANALYSIS 



Starting with the equations presented in Chapter 3, a qualitative analysis has been 
done to determine the expected trends in the numerical analysis The problem is simplified 
by assuming that the equations are exact and that there is only one satellite. The effects of 
time errors are not analyzed, so the problem is to determine how satellite position and 
focal plane error effects should effect the results 

TBM geocentric longitude and latitude can be expressed in (XYZ) coordinates: 



These (XYZ) coordinates can be expressed in terms of satellite position and focal plane 
coordinates using the LOS projection equation (8): 



Rewriting the terms on the left-hand side in terms of R, gha, 5, q, and p generates a rather 
lengthy equation, and is given in Appendix B The next task is to find the partial 
derivatives of the TBM position with respect to the five error terms, using the Chain Rule. 
Taking the partial with respect to R, for example. 



Again, expression this equation written explicitly in terms of the error sources becomes 
tedious, and can be found in Appendix B. The results, however, are interesting and are 
shown as plots. The horizontal axis pointing left is q, the horizontal axis pointing right is 
P, and the vertical axis is the differential 



X = tan — I 
Vx^ 




( 11 ) 




( 12 ) 



r= y^ =R,+pe 

UJ 



( 8 ) 



cR ox ^ dy dK dz dR 






(123) 



dz 

dR dK dR dy dR dz dR 



(124) 
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( 125 ) 



The partial derivatives with respect to gha turn out to be simple: 



5gha 




(126) 



5gha 



The observations one can make from analyzing these plots is that most of the 



partial derivative plots are flat and close to zero. When riapproaches 8.6°, however, some 
of the values begin to get large. This happens when the observed IR event is close to the 
limb of the earth as viewed from a geosynchronous satellite. There also appear to be 
transitions when r\ approaches 8 6° and b is 0° or 180°. 

Figures 44 and 45 show that satellite radius error effects are of such small 
magnitude that they are neglible Figures 46 and 47 show that satellite declination error 
effects are more important to determine latitude than longitude, except when ri 
approaches 8.6° Figures 48 and 49 show that P LOS errors affect latitude determination 
more than longitude, except, again, when ti approaches 8.6°. Figures 50 and 51 show that 
T| LOS errors have the largest effect of all, and the partials exhibit the same behavior near 
the limb. 

These plots of formulae give an indication of what magnitude the error ellipse 
should have, given the magnitude of the error is known, by using that magnitude in 
equations similar to (127); 



where R can be any of the five error quantites (R, 5, gha, r\, or P) 

To test this assumption, substitute the magnitude of the errors used in the 
simulation, and compare the result with the numerical analysis. The errors are: 




(127) 



AR = 200m 



(128) 




( 129 ) 



70 




4 I 





Agha = tan 



0200 km ^ 
.42164.17 kmj 



(130) 



Ati = 5 liradians (131) 

Ap = 5 [iradians (132) 

The generalized equations for the partial derivatives of the geocentric longitude 
and latitude can be used to find the derivatives for a specific location by using specific 
focal plane coordinates. For Baghdad, viewed from the satellite at 70° gha, in particular, 
these coordinates would be: 

11 = 0.1216 radians (133) 

P = 3.8546 radians ( 1 34) 

Using these values, the partial derivatives with respect to the five error sources 
become 



5p rad 




(135) 


^ = -0.753 
5P rad 




(136) 


^=-.3.788^ 
^ rad 




(137) 


^ = 5,96^ 
rad 




(138) 


^ ^ , rad 

— = 0.656 

rad 




(139) 


^ = -0.658^1 
55 rad 




(140) 


— = -3.663 xl Q- 
cR 


s rad 
m 


(141) 


^= 1.583 xl0 ‘ 
5R 


rad 

m 


(142) 


5X _ j rad 
5gha rad 




(143) 
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_^ = o — 

5gha rad 




(144) 


Multiplying (135) through (144) by the appropriate error quantities produces the 


geocentric longitude and latitude error values 






=3.166x10^ rad 




(145) 


A<j>;h. = 0 rad 




(146) 


A(|>; =-3.121x10^ rad 




(147) 


A(|>; =2.98x10'’ rad 




(148) 


A(|)' =-3.765x 10^ rad 




(149) 


AXr =-7.326x 10“* rad 




(150) 


=4.743x 10-^ rad 




(151) 


Ak, =-3.112x10“* rad 




(152) 


AX, = -6.894 X 10 ’ rad 




(153) 


AXp =-2 865 x 10^ rad 




(154) 


Grouping these in terms of satellite position and focal plane and summing the absolute 


values: 






A* »,o, = + A*; = 6.287 X 10-* 


rad 


(155) 


= AX, + AX,,, + AX, = 1,518 X lO’’ 


rad 


(156) 


= A*; + A*; = 3.357 x 10-' rad 




(157) 


AX„».,,u.,.„„ = AX, + AX, = 7,181 X 10 ' rad 




(158) 



These values are in radians, which can be equated to kilometers on the surface of 
the earth by multiplying latitude by the earth’s radius and longitude by the earth’s radius 
and the cosine of the latitude For this analysis, earth’s radius is assumed to be 6378.137 
kilometers: 

=(6.287x 10^ rad) 6378.137 km = 4.010x10'' km (159) 

= (1518 X lo ’ rad) cos(3 3.333“ )6378 1 3 7 km = 8.089 x 10'' km (160) 
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A(|);«.,pu.c.a« =(3.357 x 10 rad) 6378.137 km = 2.141x10-’ km (161) 

A>-roc»ipi«.ecn« = (7181 X 10-' rad) cos(33.333°)6378.137km = 3.827 x 10 ’ km (162) 
Assuming that these values are equivalent to the semi-axes of an error ellipse, multiplying 
longitude error by latitude error and by n gives the area of the error ellipses: 

ellipse area =1.019 x 10‘* km^ (163) 

ellipse area = 2.574 x lO ’ km^ (164) 

The focal plane error effects are 25 times greater than the effects of the satellite position 
error. The relative magnitudes agree with the numerical results, but the specific values 
obtained here are greater than those obtained numerically by a factor of about three The 
numerically obtained values were 

launch ellipse mean area„^„^p„i^„„ = 5.27 x 10"' km^ (165) 

launch ellipse mean areaf^.,p^„,„„ = 8.45 x 10'^ km^ (166) 

This difference is expected, since this analysis does not account for stereo viewing, the 
least squares iteration process, or time effects. 

This result shows that it is possible to compute error ellipse areas using a 
simplified qualitative analysis. The relative order of the magnitude of the effects is in 
agreement with predictions, and verifies the results obtained numerically To do the same 
analysis for burnout volumes and impact areas would require the addition of time errors, 
and trajectory equations. This increases the difficulty of this “back of the envelope” 
analysis, and is not treated here 
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VII. CONCLUSION 



The TALON SHIELD / ALERT state vector estimation algorithm vv^as correctly 
modeled in MATLAB™. System errors were modeled and introduced into the algorithm 
to determine the effects upon the accuracy of the final results. A simplified qualitative 
analysis verified the quantitative results. 

The three questions posed in the first chapter are now answered: 

1 . The errors in the system are broadly categorized as errors in time, satellite 
position, and LOS 

2. The relative magnitudes of the error effects listed largest to least: 

A LOS - using a zero-mean, 5 pradian standard deviation normal error 
distribution in focal plane coordinates 

B. satellite position - using a zero-mean, 200 meter standard deviation normal 
error distribution in satellite positon coordinates 

C. time - using a uniform error distribution between zero and 1 millisecond 

3 The effects seem to behave independently, since the principle of superposition 
works. The effects of the combined errors are slightly greater than the sum of 
the individual effects 

Additional runs were made to determine the effects of various geometries, spin rates, and 
observation synchronization 

There are many tasks left to accomplish in this area of research. The MATLAB 
code can be optimized, additional parameters and system changes can be simulated, and 
more detailed analysis of the results can be made The error sources can be modeled more 
exactly by analytically “peeling the error onion”, or by using actual DSP values and 
measurements The results obtained could be further analyzed, possibly to 

Simulating the DSP system with computer code has proven to be an effective tool 
for error analysis, which is a necessary first step for determining how best to improve the 
present TBM detection system or modify the design of a follow-on system. 
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APPENDIX A. “A” MATRIX DERIVATION 



The A matrix is made up of partial derivatives of focal plane coordinates with 
respect to the tactical parameters These elements can be determined qualitatively. Exact 
equations are not required for the iteration process to converge, so small-angle 
assumptions are made to simplify the resulting equations The benefit of using these 
approximations for the matrix elements is the ease of computation and thus a reduction of 
computing time. 

The A matrix is divided into three components, Ai, Aj, and Aj The elements of 
Ai, and A 2 will be derived in order, using MATHCAD 5 0 Plus The Aj matrix is derived 
in Chapter 3 The elements of the first column of Ai are partials of latitude, ^ : 

♦ acos cos(0) sin q -> sin(6) cos ^ q cos,a q ^ 

Taking the partial derivatives. 

jp sin(0) sin 0 q ^ cos(0) cos 1 ^ q cos a q 

d0~ “ 2 

1 cos(0) sin 0 Q - sin(0) cos 0 q cos a q 

Using the small angle assumption that sin(0)=O and cos(0)=l : 
d^ cos(|io cos ao 



2 



d a 



T T 



0 " ^2 



1 sin, (j» Q 

6 .“ 

■■efi 

d9_ I 

M rj.fl 

d=(l - 1 5L) dp 
2 



T Tn 



T Tn 



1 T Tn 
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-^ = (1 1.5 L) i - a I - 2 32 iT - Tq - 3 B 3 T - 4-34 T - Tq ^ 

dT 0 

dd ' 214 

1.5 aotaj iT- To -32.1 - Tq - 33 - T Tq -34 T-Tq 
qL 

d(f d^ d0 dd dd 

dT Q d6 dd dT Q r dT q 

^_d4 d0 dd_ '^ ^ dp cos a q 

dL" d0 dd dr’ r ^g- 

d* ^1 

dp Q 
d/. 0 

dho 



= G sin a q 
dao ^ 

The second column is composed of partials of longitude, X: 
sin(0) sin a q 
cos(<(» ) 

d'K 1 sin a Q sin a 

d0~ 



K-/. Q - 3Sm 



cos(0) 

2 COS((l>) C0S((t>) 



1 sin(G) 



'.j cos((D) 

d/. d/, d0 dd sinao- dd sin aQ. dd 

— — . . s_. — . -S ^ ^ 

dT Q d0 dd dT Q r gg'C 0 s(p ) dT q r gff cos p q dT q 
A ssuming (j> is approximately (j>o; 
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1 .5 d p sin a q 



d/. 

0 



d/._d/v d0 dd_ 1 Sdpsinao _ 

dL d9 dd dL rg^cos(O) rgg-cos> q 
1 



do 

sin(G) - - -^ sin(4) =sin(0) sin a q 



1 sin(0)^ 



‘ o ; 



cos(^) 



cos(o) 



d?. _ 

d ?. Q 

dhr 



=0 



d *( 



da 



0 



- sin(0) = - 

2 cos((+)) cos^^ Q 



1 sin(G) 



cos(^r 

And the third column elements are the partials of height, h; 

h*ho (1 - L) hp 



'’p'l’O *>1 T To-b2T Tg^bjT To^-b4.T 

dh f 2 

s(l + L) b , 2 b 2' 7 Tq - 3by J - 7q - 4 b ^ T 

dTo 



‘'^h 

dL 

dh _ 

dp Q 

dh _ 
d/. n 



dh _ 
dho "^ 



sin(o) 

cos(p 
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dh 
da 0 

dh^ 

dho' 



The A 2 matrix is the partials of (XYZ) with respect to (<J>, X, h): 
Xt, = cos(4)\ ) + altk cos(<|>j]cos(X J 

yj, = cos( ({)',) + alt k cos((J>j]sin(X,) 

= rfoci, sin(4>\ ) + alt^ sin(<|>J 

Approximating riocai with rcff , alt with h, (}>' with and taking the partials: 
-^ = -(r^ ^h)sin(4))cos(X) 

-^ = -(fcff +h)cos((j>)sin(X) 



— = cos(4>)cos(X) 
(Ih 



-^ = -(r^ -h)sin(<f))sin(X) 
ext) 



— = (r^g + h)cos((j))cos(X) 
c;X 



= cos(4))sin(X) 
dh 



dz 



(r,jj +h)cos(4>) 



^=0 

a 



— = sin((j>) 
ai ^ 
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APPENDIX B. QUALITATIVE ANALYSIS EQUATIONS AND PLOTS 
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o. i c- 
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iR cos(5) - rioj-ai 
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The geocentric longitude and latitude can be expressed in terms of (XYZ) 
coordinates 



lx' 

^ 



1 2: 

X -1 

1 1 


u>Ll| 

x/j 

1 


1 ' 


-Si 

z ' 

1^1 

X 


dx ,3 






1 + ! 




y 


dy 


: 




I 2 2 

1. X T >■ 


dU_ 


1 



Vx - y • 1- 



The differentials of the geocentric longitude and latitude can be determined using 
the Chain Rule Using longitude and satellite radius as example quantities: 
dX _ dX dx dX dy 

Substituting the following quantities into the previous equations: 

R =42164 17kn) 
gha =70 deg 
8 =0 

r local =6378.137kn) 

Ti = vanes O' to S.S' 
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P = varies 0° to 360° 



This gives expressions for the variations of the geocentric longitude and latitude 
due to errors in LOS projections from a DSP satellite with a Greenwch hour angle of 70°. 
These expressions are transformed into plots for examination. 

The following plots show the differentials of geocentric longitude and latitude with 
respect to the five error sources. The horizontal axis pointing to the left is r\, the 
horizontal axis pointing toward the right is 3, and the vertical axis is the differential 
evaluated at each (t|, P) coordinate. 
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d(j)'dp. ^ - dijiUx. j dxdp. + d(t>'dy. ydydp, . + dfdz. j dzdp. ^ 




d^'dp 




dXdp 
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dfdti. j = d(J'd.\. j dxdtij ^ -t- dfdy. j-dydr]. ^ + d^'dz. j dzdti. ^ 





d>.dii 
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dA.d5. . =dA.dx. . dxd5 . + dA.dy . dydS. . 

i,j i,j i,j ^ i,j 




dXd5 
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dfdR, j = dfdx. j dxdR. . + d(|('dy. j dydR. ^ + diji'dz. /dzdR. . 




d(|i'dR 



dldR. . = dXdx. y dxdR. . + cEdy. . dydR. . 




d^LdR 
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tlijidgha. j - diMx. ydxdgha. ^ +- d(tidy. . dydgha. ^ 




d(^gha 



dXdgha. j =<D.dx. ydxdglia. . 4- d).dy. y dydgha. ^ 




d^gha 
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nJA Aj?'- 
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APPENDIX C. MATLABCODE 



The following is the MATLAB™ program used to simulate the TALON SHIELD 
/ ALERT state vector estimation algorithm. It is commented in blue font. 

®/o Martin Beaulieu 

% DSP Estimation .Algorithm 



9o Initialize variables 




clear; 


% clear memor\ 


format long, 


% use 1 5-digit numerical format 


numpoints=1000. 


% number of data points to run 


reff^6371, 


% spherical eanh radius (km) 


f=l/298.257. 


°b WGS-84 flattening factor for oblate earth 


re=6378.137; 


°-o oblate earth equatorial radius (km) 


a0=0.074929 1380897402, 


% coefficients in range profile equation 


a 1 =-0.050 164 7424642791; 




a2=0.004 194795 19493 87, 




a3=-0.0000141392516462238. 




a4=8.51293969732759E-07; 




b0=0.884282320240989. 


“o coefficients in height profile equation 


bl=-0.129163814041924. 




b2=0.01 37288843548327, 




b3=-0.0001 59307288768991, 




b4= 1 .36938557528 199E-06, 




tmax=62.5. 


°o TBM maximum motor burn time (second? 


mu=3 9860 1.2, 


°o earth gravitational constant (km '3 sec '2) 


d2r=pi/180. 


°o degrees to radians 


r2d= 180/pi; 


°o radians to degrees 


we=15*d2r/3600, 

0 


°o earth rotation rate ( 15 degrees 'hour) 



®'o This section generates the observ ation table (time. az. el. gha. dec. and radius) The first 
°o step is to choose the tactical parameters, then calculate hovs that missile profile would 
% appear to an orbiting sensor. 

% Initial tactical parameters 

for ii=0;l 1 °/o loop through 12 headings around compass rose 

T0=100; ° o launch time of day in seconds 

L=0, % loft parameter 

%latO=l 2.783 *d2r;lon0=45.050*d2r, % Aden, A emen 
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Iat0=3 3 . 3 3 3 * d2r,lon0=44 .43 3 * d2r, 
%latO=3 3 . 500 * d2r;lonO=3 6. 3 1 9 * d2r, 
%lat0=24.633 *d2r;lon0=46. 7 1 7*d2r; 
%lat0=3 5.666*d2r;lon0=5 1 .433 *d2r; 
h0=0; 

hdg0=ii*pi/6, 

data=[T0 L lonO latO hO hdgO zeros(l,9)], 

% Generate obser\ation time matrix- 

obsl =T0+30+5*randn( 1 ), 
dtl=10*rand(l),dt2=10*rand(l), 
ttime=[obs 1 ;obs 1 +dt 1 ;obs 1 +dt2]; 
ttime=sort(ttime), 
nextobs=obsl+10; 

j=n 

while nextobs <= tmax+TO, 
ttime=[ttime,ttime(j)+ 1 0], 

nextobs=ttime(j )+ 1 0 , 
end 

n=length(ttime), 

“o (ienerate satellite position matrix 

gha=[ 1 0*d2r,70*d2r,l 05*d2r], 
satord=ceil(6*rand(l )); 
if satord = 1 

tgha=[gha(l),gha(2),gha(3)], 
elseif satord = 2 
tgha=[gha( 1 );gha(3),gha(2)], 
elseif satord = 3 
tgha=[gha(2),gha( 1 ),gha(3)], 
elseif satord = 4 
tgha=[gha(2);gha(3),gha( 1 )]; 
elseif satord = 5 
tgha=[gha(3),gha( 1 ),gha(2)], 
elseif satord = 6 
tgha=[gha(3);gha(2),gha( 1 )], 
end 

tdec=[0,0;0]; 

tradius=[42 1 64. 1 7,42 1 64. 1 7,42 1 64. 1 7], 
for i=4:n 

tgha=[tgha,tgha(i-3)]; 

tdec=[tdec;tdec(i-3)]. 



% Baghdad. Iraq 
% Damascus, S\Tia 
% Riyadh. Saudi .Arabia 
% Tehran. Iran 
% launch altitude 

% launch heading in 30 degree steps 
% true tactical parameters 



° o first observ ation time - clouds, vapor, etc 
% two random time intervals. 0<dt<10 sec 
% first 3 observation times 
% put times in chronological order 
% next sequential observation 
counting index 

% if next observation is during motor burn 
add next obseix ation to time matrix 
°o increment counter 
°o next sequential observ ation 
°o loop until matrix is completed 
°o number of observations 



“n (ireenwich hour angle 
°o randomlv select the order in 
°o which .satellites observe TBM 



°o declination 

°o geosv nchronoLis radius 

% make satpos matrix the same size as time 
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tradius=[tradius;tradius(i-3,:)], 

end 

satpos=[tradius. *cos(tdec). *cos(tgha) 
tradius *cos(tdec). *sin(tgha) 

tradius *sin(tdec)]; °/o satellite position 

% Apply profile to observation times: 



t=ttime-T0, 

dp=aO+a 1 *t+a2 *t /^2+a3 *t /'3+a4 *t . 
d=(l-1.5*L)*dp; 

hp=b(H-bl *t+b2*t.^2+b3*t.^3+b4*t 
h=(l+L)*hp; 



% time of flight 
% nominal profile distance 
% loft-modified distance 
nominal profile height 
% loft-modified height 



% Generate target position: 



theta=d/reff, % eanh-central angle 

lati=pi/2-acos(cos(theta)*sin(latO)+sin(theta)*cos(latO)*cos(hdgO));® o latitude 
long=lonO+asin(sin(theta)*sin(hdgO)./cos(lati)); % longitude 

alt=hO+h, % altitude 

gcl=atan((l-f)'^2*tan(lati)), °-o geocentric latitude 

rloc=re*(l-f)./sqrt((l-f)^2*cos(gcl)/'2+sin(gcl)/'2), ®o local radius 

tgtpos=[(rloc.*cos(gcl)-t-alt.*cos(lati)).*cos(long) 

(rloc. *cos(gcl)+alt. *cos(lati)). *sin(long) 

rloc.*sin(gcl)+alt.*sin(lati)], % target position 



°o Generate line-of-sight sector and transform into focal plane coordinates 



delta=tgtpos-satpos, °-o line-of-sight vector 

r3=[];r3t=[],UEN=[]; ® o reset variables 

for i=l :n 

rotl=[cos(tgha(i)) -sin(tgha(i)) 0,sin(tgha(i)) cos(tgha(i)) 0;0 0 1]; 

rot2=[cos(tdec(i)) 0 -sin(tdec(i));0 1 0;sin(tdec(i)) 0 cos(tdec(i))]; 

rot3=rotl *rot2; 

r3t=[r3t,rot3'], 

r3=[r3;rot3]; 

UEN(i,:)=(r3t(3*i-2:3*i,:)*delta(i,:)')', % Up-East-North coordinates 

end 

tbeta=rem(-pi/2-atan2(UEN(:,3),UEN(:,2))+(2*pi),(2*pi)),°o true azimuth 
teta=atan(sqrt(UEN(:,2)/'2+UEN(:,3)/^2)./(-UEN(;,l))); true elevation 



®o L sing generated data (ttime, theta, teta. tgha. tdec. tradius), 

®o compute tactical parameters using least-squares iteration method 
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for index= 1 : numpoints, 
index=round(index); 



% ensure index is an integer 



% Initialize variables 



clear time radius gha dec beta eta TO L latO lonO hO altO hdgO tp; 
sse=l ;dsse=l ;ssel=l ,failindex=0,r3t=[];r3=[]; % reset variables 

errtime=zeros(n, 1 ), % time error 

enrad=zeros(n, 1 ); “/o satellite radius error 

errgha=zeros(n,l), % satellite gha error 

errdec=zeros(n, 1 ); % satellite dec error 

errbeta=zeros(n, 1 ); % azimuth error 

erreta=zeros(n,l); % elevation error 



°-o Change time, beta. eta. and satpos by some error, epsilon 



if index > 1 

errtime= 1 e-3 *(rand(n, 1 )-0. 5), 
errrad=0.2*randn(3,l ); 
errgha=atan(0.2/42 1 64. 1 7)*randn(3, 1 ), 
erTdec=atan(0.2/42164. 1 7)*randn(3,l), 
for i=4;n 

enTad=[errrad,errrad(i-3 )] , 
errgha=[eiTgha;errgha(i-3 )], 
errdec=[errdec ,errdec(i-3 )] , 
end 

errbeta=5e-6*randn(n, 1 ), 
erreta=5e-6*randn(n, 1 ), 
end 

time=ttime+errtime, 

radius=tradius+errrad, 

gha=tgha+errgha, 

dec=tdec+errdec; 

beta=tbeta+errbeta, 

eta=teta+erreta; 



°o <=--0 Oul second 
® 0 variance of -^-200 meters 
°o variance of -^-200 meters 
°o variance of *-200 meters 
°o fix satellite positions 

“ o during single run 



v ariance^ --Smicroradians 
°o \ ariance= - -Smicroradians 

°o time = truth - error 
°o radius = truth - error 
gha = truth * error 
% dec = truth * error 
°o beta = truth error 
eta ^ tiTJth - error 



% Compute target position. 

satpos=[radius. *cos(dec). *cos(gha) 
radius. *cos(dec).*sin( gha) 

radius. *sin(dec)]; % satellite position 

for i=l :n % rotation matrices 

rotl=[cos(gha(i)) -sin(gha(i)) 0;sin(gha(i)) cos(gha(i)) 0;0 0 1]; 
rot2=[cos(dec(i)) 0 -sin(dec(i));0 1 0;sin(dec(i)) 0 cos(dec(i))]; 
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rot3=rotl *rot2, 
r3t=[r3t;rot3'], 
r3=[r3;rot3]; 
uen=[-cos(eta(i)), 

-sin(beta(i )) * sin(eta(i )), 

-cos(beta(i))* sin(eta(i))] , 
ehat=rot3*uen, 

chi=pi-asin((norm(satpos(i,:))*sin(eta)/reflO), 
rho=reS7sin(eta) * sin(eta+chi); 
los=rho*ehat; 
tgtpos=los+satpos(i , ; 
lat(i)=atan(tgtpos(3 )/sqrt(tgtpos( 1 )'^2+tgtpos(2)'^2)); 
lon(i)=atan2(tgtpos(2),tgtpos( 1 )), 
end 



% Up-Easi -North coordinates 
% line of siglit direction 
°'o computation of length 
°/o for los vector 
line-of-sight vector 
° o target position 
°'o latitude 
°-b loncitude 



°'o Intial estimate of tactical parameters 

lat0==(lat(l)+lat(2)+lat(3))/3, °o first obs latitude 

lon0=<lon(l)+lon(2)+lon(3))/3, ° o first obs longitude 

latbo=(lat(n-2)+lat(n-l)+lat(n))/3, °o latitude at last obs 

lonbo=(lon(n-2)-t-lon(n-l)+lon(n))/3, °o longitude at last obs 

hdg0=rem(pi/2-atan2(latbo-lat0,(lonbo-lon0)*cos(lat0))+(2*pi),(2*pi)), ° olnch heading 
T0=time(l)-20; launch time 

L=0, ° 0 loft pai ameter 

h0=0; ° o launch height abo\c WGS-84 ellipsoid 

tp=[T0^,lat0don0,h0^dg0], °o tactical parameters matri\ 

% Iterate on initial estimates until the difference between the sum of the 

% squares of the errors between two consecutive runs is <=10*eps: 



while dsse > 10*eps 
duv=[],A=[], 

T0=tp(l );L=tp(2),latO=tp(3); 

Ion0=tp(4),h0=tp(5),hdg0=rem(tp(6),2 * pi), 
t=time-T0; 

dp=aO+al *t+a2*t.^2+a3*t."'3+a4*t.^4, 
d=(l-1.5*L)*dp, 

hp=bO+b 1 *t+b2 *t.'^2+b3 *t ,^3+b4 *t.M, 
h=(l+L)*hp, 
theta=d/reff, 

lati=pi/2-acos(cos(theta)*sin(latO)+sin(theta)*cos(latO)*cos(hdgO));°o geodtic latitude 
long=lonO+asin(sin(theta)*sin(hdgO)./cos(lati)), % longitude 
alt=hO+h, °'o altitude 

gcl=atan((l-f)'^2*tan(lati)), °o geocentric latitude 

rloc=re*(l-f) /sqrt((l-f)^2*cos(gcl).^2+sm(gcl)/^2);°b local radius 



°o loop until dsse '' 10*eps 
°o reset \anables 
° o update tp matrix \ alues 
° o with dtp's added 
°'o time-of-flight 
°o nominal profile distance 
°/b loft-modified distance 
% nominal profile height 
% loft-modified height 
°'o eanh-central angle 



97 



tgtpos=[(rloc. *cos(gcl)+alt. *cos(lati)). *cos(long) 
(rloc. *cos(gcl)+alt. *cos(lati)) *sin(long) 
rloc.*sin(gcl)+alt.*sin(lati)]; ° 



*?'0 target position 
°o line-of-sight vector 
% rotate into UEN coordinates 



delta=tgtpos-satpos, 
for i=l:n 

UEN(i,:)=(r3t(3*i-2;3*i,:)*delta(i,:)')', 

end 

uhat=-UEN(;,2)./UEN(:,l), 

vhat=-UEN(;,3)./UEN(:,l); 

u=-tan(eta). * sin(beta), 

v=-tan(eta). *cos(beta); 

du=u-uhat, 

dv=v-vhat. 



°o estimated focal 
% plane coordinates 



°'o difference bet \^een 

% estimate and actual 



actual focal 
plane coordinates 



°'o Compute A matrix: 

cl=cos(hdgO)/reff,c2=sin(hdgO)/refi7cos(latO); constants in A matrix 

c3=sin(hdgO)/reff,c4=cos(hdgO)/reff7cos(latO), 

dddT0=-(l-1.5*L)*(al+2*a2*t+3*a3*t.''2+4*a4*t.'^3), °o -horizontal speed 
dhdT0=-(l+L)*(bl+2*b2*t+3*b3*t.^2+4*b4*t/'3); °o -vertical speed 
dudUEN=[UEN(;,2)./UEN(:,l).^2 -l./UEN(;,l) zeros(size(time))]; 
dvdUEN=[UEN(:,3)AJEN(:,l)/'2 zeros(size(time)) -l./UEN(:,l)], 

°o change in local plane coordinates \ui changes in I EN coordinate^> 
for i=l;n construct A matrix 

duv=[duv;du(i),dv(i)], " c error \ allies matrix 

Al= [dddT0(i)*cl dddT0(i)*c2 dhdTO(i), 

-1.5*dp(i)*cl -1.5*dp(i)*c2 hp(i), 

1 0 0,0 1 0;0 0 l;-d(i)*c3 d(i)*c4 0]; 
A2=[-(rloc(i)+alt(i))*sin(lati(i))*cos(long(i)) 
-(rloc(i)+alt(i))*sin(lati(i))*sin(long(i)) 

(rloc(i)+alt(i))*cos(lati(i)); 

-(rloc(i)+alt(i))*cos(lati(i))*sin(long(i)) 
(rloc(i)+alt(i))*cos(lati(i))*cos(long(i)) 0; 
cos(lati(i))*cos(long(i)) cos(lati(i))*sin(long(i)) sin(lati(i))], 
Atmp=Al*A2*r3(3*i-2:3*i,:); % A 3 -(delta ' delta x>z)-r 

A=[A;(Atmp*dudUEN(i,;)')',(Atmp*dvdUEN(i,;)')'], ° o total A matrix 



end 



AT=pinv(A); 

dtp=AT*duv; 

tp=tp+dtp; 

ssel=sum(duv/'2); 

dsse=abs(sse-ssel); 

if ssel < sse 



% find pseudoinverse of A 

°o find changes to tp 

° 0 add changes to tp 

°b find sum of the squares of the error 

®-o difference between iterations 



sse=ssel; 

end 
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failindex=failindex+ 1 , 
if failindex>100 
dsse=0, 
end 
end 

% Burnout lime estimation 
Tlast=time(n), 

Tnext=2*time(n-2)-time(n-5), 
Tmax=TO+tmax, 
if index =1 
Tbo=T0+tmax, 
el seif T max>=Tnext 
Tbo=Tlast-K) . 5 * (T next-Tlast), 
else 

Tbo=Tlast +0 . 5 * (T max-Tlast ), 
end 

tbo=Tbo-T0, 

® o Calculate^ state \ector at burnout 



% increment counter 

% if statement to prevent continuous looping 
% in the event that it does not converge 



°/o last obsersation time 

°o next obser\ ation if TBM was still burning 

% max observation time from profile 

° o pick bum out time as half of the time 
°/o between Hast and Tnext or Hast and 
° 'o Tmax depending on relative 
°o magnitude of Tmax and Tnext 

°o burnout time-of-flight (close to tmax) 



% loop to next tp iteration 



d=(l-l 5*L)*(a0+al*tbo+a2*tbo''2+a3*tbo''3+a4*tbo''4), “o distance 

ddot=(l-1.5*L)*(al+2*a2*tbo+3*a3*tbo^2+4*a4*tbo^3), horizontal speed 

h=(l+L)*(b0+bl*tbo+b2*tbo^2+b3*tbo^3+b4*tboM), height 

hdot=(l+L)*(bl+2*b2*tbo+3*b3*tbo'^2+4*b4*tbo^3), vertical speed 

theta=d/reff, ® o eanh-central angle 

latbo=pi/2-acos(cos(theta)*sin(latO)+sin(theta)*cos(latO)*cos(hdgO));® ogeodtic latitude 
lonbo=lonO+asin(sin(theta)*sin(hdgO)/cos(latbo)), ® o burnout longitude 

hbo=h0+h, % burnout altitude 

Vbo=sqrt(ddot'^2+hdot'^2); ® o burnout velocitv 

fpabo=atan(hdot/ddot), ® o flight path angle 

hdgbo=rem(asin(cos(lat0)*sin(hdg0)/cos(latbo))+(2*pi),(2*pi)), ° o burnout heading 



°o Calculates ballistic trajectorv of target to impact 

altbo=h0+hbo, % burnout altitude 

gclbo=atan(((l-f)^2)*tan(latbo)), geocentric latitude 

rlocbo=Te*(l-f)/sqrt((l-f)^2*cos(gclbo)^2+sin(gclbo)^2), % local radius 

tgtposbo=[(rlocbo*cos(gclbo)+altbo*cos(latbo))*cos(lonbo) 

(rlocbo * cos(gclbo)+altbo * cos(latbo))* sin(lonbo) 
rlocbo*sin(gclbo)+altbo*sin(latbo)]; % TBM position at burnout 
rbo=norm(tgtposbo), % radius at burnout 

vo=Tlocbo*we*cos(gclbo), ® o initial inertial velocity 

Vu=Vbo*sin(fpabo), “b up component of v eiocitv 
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% east component of velocity 
% north component of velocity 
% inenial speed 
% inertial flight path angle 
% inertial heading 
°/o energy' parameter 
% eccentricity of orbit 
% freeflight eca 
% eccentric anomaly 



Ve=Vbo*cos(fpabo)*sin(hdgbo)+vo, 

Vn=Vbo*cos(f^abo)*cos(hdgbo), 

Vin=sqrt(Vu/^2+Ve.^2+Vn.^2); 
fpain=asin(VuA^in), 

hdgin=Tem(pi/2-atan2(-Ve,Vin)+2*pi,2*pi); 

Qin= Vin'^2 * rbo/mu; 
em=sqrt( 1 +Qin*(Qin-2)*cos(fpain)^2); 

PSI=2*acos(( 1 -(^*cos(fpain)^2)/ein)+theta, 

E=acos((ein-cos(PSI/2))/( 1 -ein*cos(PSI/2))); 
ain=rbo/(2-Qin), % semimajor axis 

tfi=2*sqrt(ain'^3/mu)*(pi-E+ein*sin(E)); % time of free flight 

gclim=asin(sin(gclbo)*cos(PSI)+cos(gclbo)*sin(PSI)*cos(hdgin)); %geocentric latitude 
latim=atan(tan(gclim)/(l-f)^2); % geodetic latitude of impact 

Cnon=acos((cos(PSI)-sin(gclim)*sin(gclbo))/(cos(gclim)*cos(gclbo)))-we*tff, 
lonim=lonbo+fflon, % longitude of impact 

tim=2*tbo+tff, °'o total flight time 

Tim=T0+tim, °/o time of day of impact 

data=[data,tp(l) tp(2) tp(4) tp(3) tp(5) tp(6) lonbo latbo hbo fpabo hdgbo Vbo Tim 



lonim latim], 
end 



°o sa\e data in matrix 
°'o end of numpoinl loop 



°o Plotting routines and data anahsis 

tllon=data(2,3 ) * r2d; 
tllat=data(2,4)*r2d, 
llon=data(3 . index+ 1 ,3 )* r2d, 
llat=data(3 : index+ 1 ,4)*r2d, 
mllon=mean(data(3 : index+ 1 ,3 ))*r2d, 
mllat=mean(data(3 :index+ 1 ,4))*r2d, 
dllon=tlJon-mllon; 
dllat=tUat-mllat; 



®o true launch longitude 
°o true launch latitude 
“o corrupt launch longitude 
° 0 corrupt launch latitude 
° o mean of corrupt longitudes 
°/o mean of corrupt latitudes 
°o diff between mean and true lion 
°o dilT between mean and true Hat 



tblon=data(2,7)*r2d, 
tblat=data(2,8)*r2d; 
tbalt=data(2,9); 
blon=data(3 ; index+ 1 ,7)*r2d, 
blat=data(3 : index+ 1 ,8)* r2d, 
balt=data(3 ;index+ 1 ,9); 
mblon=mean(data(3 : index+ 1 ,7))*r2d; 
mblat=mean(data(3 : index+ 1 ,8 ))* r2d, 
mbalt=mean(data(3 : index+ 1 ,9)); 
dblon=tblon-mblon; 
dblat=tblat-mblat, 
dbalt=tbalt-mbalt, 



°/o true burnout longitude 
° o true burnout latitude 
% true burnout altitude 
°'o corrupt burnout longitude 
% corrupt burnout latitude 
% corrupt burnout altitude 
% mean of corrupt bo longitudes 
°/'o mean of corrupt bo latitudes 
% mean of corrupt bo altitudes 
diff betw een mean and true blon 
°o dilT between mean and true blat 
° o diff betw een mean and true blat 
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tilon=data(2, 14)*r2d; 
tilat=data(2,15)*r2d, 
ilon=data(3 : index+ 1 , 1 4) * r2d, 
ilat=data(3: index+ 1 , 1 5)*r2d, 
milon=mean(data(3 : index+1 , 14))*r2d, 
miJat=mean(data(3 ; index+ 1 , 1 5)) *r2d; 
dilon=tilon-milon, 
dilat=tilat-milat, 

lpts=[llon-tllon llat-tllat]; 
lk=cov(lpts), 

[lv,ld]=eig(lk), 
lab=2 * sqrt(diag(ld)), 
la=lab(l); 
lb=lab(2), 

lth=pi/2-atan(lv( 1 , 1 )/lv(2, 1 )), 

bpts=[blon-tblon blat-tblat balt-tbalt], 
[bv 1 ,bd 1 ]=eig(cov(bpts(: , 1 ),bpts(: ,2))), 
babc=2*sqrt(diag(bd 1 )), 
ba=babc( 1 ),bb=babc(2),bc=babc(3); 

ipts=[ilon-tilon ilat-tilat], 
ik=cov(ipts), 

[iv,id]=eig(ik); 
iab=2 *sqrt(diag(id)), 
ia=iab(l); 
ib=iab(2), 

ith=pi/2-atan(iv( 1 , 1 )/iv(2, 1 )), 

°o Sa\ e data to anahze 



% true impact longitude 
% true impact latitude 
% corrupt impact longitude 
% corrupt impact latitude 
% mean of corrupt im longitudes 
% mean of corrupt im latitudes 
% diff between mean and true ilon 
diff between mean and true ilat 

% center data about launch truth 

% covariance of scatter 

°'o eigen value & vector of covariance 

“ o dimensions of ellipse 

°o launch ellipse semimajor axis 

°o launch ellipse semiminor axis 

°o angle to rotate ellipse 

°o center data about burnout truth 
®o eigen val and vector of covariance 
dimensions of ellipse 
"o burnout ellipse semimajor axis 

°o center data about launch truth 
“o cov ariance of scatter 
°o eigen value S: vector of covarance 
°o dimensions of ellipse 
°o impact ellipse semimajor axis 
°o impact ellipse semiminor axis 
°o angle to rotate ellipse 



data=[data,la lb 1th bal bbl bthl ba2 bb2 bth2 ba3 bb3 bth3 ia ib ith, mllon mllat mblon 
mblat mbalt milon milat dllon dllat dblon dblat dbalt dilon dilat 0]; 
ifii=0 

save bOt data -ascii -double 
elseifii=l 

save b3t data -ascii -double 
el seif ii==2 

save b6t data -ascii -double 
elseif ii=3 

save b9t data -ascii -double 
elseifii=4 

save bl2t data -ascii -double 
elseif ii=5 
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save bl 5t data -ascii -double 
elseif ii==6 

save bl8t data -ascii -double 
elseif ii==7 

save b21t data -ascii -double 
elseif ii==8 

save b24t data -ascii -double 
elseif ii=9 

save bt27 data -ascii -double 
elseif ii= 10 

save bt30 data -ascii -double 
elseif ii=l 1 

save bt33 data -ascii -double 
end 

clear data 
end 
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